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1  Introduction 


The  main  objective  of  this  project  was  the  development  of  an  innovative  computer  program  with 
advanced  capabilities  for  providing  means  for  establishing  reliable  quantitative  failure  initiation  crite¬ 
ria  for  laminated  composites  and  adhesively  bonded  joints.  The  existing  software  product  Stress 
Check,  which  is  based  on  the  p-  and  hp- version  of  the  finite  element  method  and  has  a-posteriori  error 
estimation  capabilities,  provided  the  framework  for  this  development. 

The  evaluation  of  strength  was  based  on  recent  technological  advances  which  made  it  possible  to 
determine  the  natural  straining  modes  and  their  intensities  at  singular  points  associated  with  multi¬ 
material  interfaces  by  numerical  methods.  The  project  was  concerned  with  proof  of  concept  and 
implementation  in  a  fully  three-dimensional  setting,  in  conjuction  with  the  hp-version  of  the  finite 
element  method.  During  the  Phase  II  project,  a  set  of  failure  criteria  for  composite-laminated  bonded 
structures  based  on  stress  averaging  and  material  nonlinear  analysis  were  identified  as  important  fail¬ 
ure  prediction  methods  of  interest  to  the  Air  Force,  Navy  and  their  contractors.  Therefore  a  change  in 
the  original  technical  objectives  of  Phase  II  was  incorporated  during  the  last  6  months  of  the  project. 

The  project  utilized  an  existing  finite  element  software  product  called  Stress  Check  which  was  devel¬ 
oped  by  the  small  business  concern  (ESRD)  over  the  past  ten  years.  The  research  institution  (Wash¬ 
ington  University)  provided  specialized  expertise  needed  for  effecting  transfer  of  the  new  technology 
to  professional  practice.  The  new  technological  capability  developed  under  this  project  makes  it  pos¬ 
sible  to  improve  the  design  from  the  point  of  view  of  durability.  It  will  also  make  it  possible  to  prop¬ 
erly  interpret  experimental  observations  concerning  failure  initiation  events  in  composite  materials 
and  bonded  joints  and  thereby  advance  the  use  of  composite  materials  technology  in  both  civilian  and 
military  applications. 
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2  Technical  objectives 


During  the  Phase  II  project,  the  following  tasks  were  completed: 

Task  1:  Development,  implementation,  testing  and  documentation  of  an  algorithm  for  the  computa¬ 
tion  of  eigenpairs  that  characterize  the  temperature  distribution  in  the  vicinity  of  singular 
edges  in  multi-material  interface  for  steady  state  heat  transfer  problems. 

Task  2:  Development,  implementation,  testing  and  documentation  of  an  algorithm  for  the  computa¬ 
tion  of  the  generalized  flux  intensity  factors  (GFIFs)  for  edge  singularities  in  the  case  of 
steady  state  heat  transfer  problems. 

Task  3:  Testing  and  documentation  of  the  procedures  described  in  tasks  1  and  2.  These  procedures 
were  implemented  in  the  finite  element  analysis  software  product  Stress  Check  and  a  posteri¬ 
ori  error  estimation  capability  in  terms  of  the  flux  intensity  factors  is  provided. 

Task  4:  Development,  implementation,  testing  and  documentation  of  an  algorithm  for  the  computa¬ 
tion  of  eigenpairs  that  characterize  the  strain  distribution  in  the  vicinity  of  singular  edges  in 
multi-material  interface  elasto-static  problems. 

Task  5:  Development,  implementation,  testing  and  documentation  of  an  algorithm  for  the  computa¬ 
tion  of  the  generalized  stress  intensity  factors  (GSIFs)  for  edge  singularities  when  an  elastic 
body  is  subjected  to  mechanical  and  thermal  loads. 

Task  6:  Testing  and  documentation  of  the  procedures  described  in  tasks  4  and  5.  These  procedures 
were  implemented  in  the  finite  element  analysis  software  product  Stress  Check  and  a  posteri¬ 
ori  error  estimation  capability  in  terms  of  the  stress  intensity  factors  is  provided. 

Task  7:  Development,  implementation,  testing  and  documentation  of  a  capability  for  the  computation 
of  average  stress/strain  along  element  edges,  element  faces  or  arbitrary  curves  for  two-  and 
three-dimensional  elasticity  problems.  The  average  is  understood  in  the  integral  sense  over  a 
characteristic  length/area. 

Task  8:  Development,  implementation,  testing  and  documentation  of  an  automatic  nonlinear  solution 
procedure  for  the  evaluation  of  composite  bonded  joints.  The  procedure  involves  a  material 
nonlinear  analysis  of  the  joint  for  the  design  load  followed  by  the  determination  of  the  failure 
load  based  on  several  predefined  criteria.  This  automatic  procedure  is  accessible  from  within 
the  handbook  framework  of  Stress  Check. 

Task  9:  Development,  implementation,  testing  and  documentation  of  a  capability  to  input  orthotropic 
material  properties  for  individual  plies  and  for  sub-laminate  properties  for  2D  plane-strain 
analyses.  The  procedure  allows  entering  the  3D  material  coefficients  in  the  material  axes  of 
the  composite  and  performs  the  necessary  transformations  to  compute  the  equivalent  2D 
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properties  in  the  Stress  Check  XY  work  plane.  When  a  set  of  plies  need  to  be  combined  in  a 
single  layer  (sub-laminate),  the  properties  of  the  sub-laminate  are  determined  by  homogeni¬ 
zation  of  the  material  properties. 

Task  10:  Delivery  of  one  copy  of  Stress  Check  software  executable  on  a  Windows  NT  workstation 
and  pertinent  documentation  to  the  Air  Force  Rome  Laboratory.  Rome  Laboratory  shall  have 
the  same  rights  to  Stress  Check  as  a  “paid-up”  licensee. 

The  details  of  these  activities  are  described  in  the  following  sections. 

3  Technical  background 


There  is  a  growing  demand  for  efficient  and  reliable  means  for  predicting  and  eventually  prevent¬ 
ing  failure  initiation  and  propagation  in  multi-chip  modules  (MCM),  electronic  packages  and  com¬ 
posite  materials  subjected  to  mechanical  and  thermal  loads.  Thermal,  elastic  and  thermo-elastic 
problems  associated  with  large  scale  integrated  circuits,  electronic  packaging,  and  composites 
increase  in  complexity  and  importance.  These  components  are  assemblages  of  dissimilar  materials 
with  different  thermal  and  mechanical  properties.  The  mismatch  of  the  physical  properties  cause  flux 
and  stress  intensification  at  the  comers  of  interfaces  and  can  lead  to  mechanical  failures. 

The  traditional  finite  element  analysis  of  stresses  is  inadequate  to  handle  these  types  of  problems1: 
“ Since  the  stress  and  displacement  fields  near  a  bonding  edge  show  singularity  behavior,  the  adhesive 
strength  evaluation  method,  using  maximum  stresses  calculated  by  a  numerical  stress  analysis,  such 
as  the  finite  element  method,  is  generally  not  valid". 

These  material  interfaces,  as  well  as  crack  tips,  are  called  singular  points  because  the  temperature  flux 
is  infinity  in  the  linear  theory  of  steady-state  heat  conduction  and  so  are  the  stresses  in  the  linear  the¬ 
ory  of  elasticity.  Typical  singular  points  where  failures  initiate  and  propagate  in  an  electronic  device 
are  shown  in  Figure  1. 

New  approaches  to  predicting  the  initiation  and  extension  of  delamination  in  plastic  encapsulated  LSI 
(Large  Scale  Integrated  Circuit)  devices,  for  example,  are  based  on  the  computation  of  certain  func¬ 
tionals,  called  the  generalized  flux/stress  intensity  factors  (GFIFs/GSIFs),  the  strength  of  the  stress 
singularity,  and  in  thermo-elastic  problems  the  thermal  stress  intensity  factors  (TSIFs).  Since  the 
stress  and  strain  components  are  generally  not  bounded,  it  is  not  possible  to  construct  failure  initiation 
events  with  them.  A  key  requirement  for  formulating  failure  laws  is  a  software-tool  which  computes 
reliably  the  GFIFs/GSIFs  and  the  strength  of  the  singularity.  The  reliability  and  accuracy  of  numerical 
results  is  an  essential  requirement  for  proper  evaluation  and  formulation  of  failure  theories  through 
experimental  observation  of  failure  initiation  events. 


1.  T.  Hattori,  S.  Sakata  and  G  Murakami,  “A  stress  singularity  parameter  approach  for  evaluating  the  interfacial  reliability  of  plastic 
encapsulated  LSI  devices",  Jour.  Electronic  Packaging,  111,  (1989),  pp.  243-248. 
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FIGURE  1.  Typical  sites  of  failure  initiation  in  an  electronic  device. 


3.1  Introduction  and  Notation 

In  the  neighborhood  of  singular  points  the  exact  solution  of  two-dimensional  elastostatic  prob¬ 
lems  for  example,  can  be  expanded  in  the  form: 


{«}  =  £  a/* {*.&)}  (D 

i  —  1 

where  {«}  is  the  displacement  vector  with  components  ux(x,y),  uy(x,y),  r  and  6  are  polar  coordinates 
centered  on  the  singular  point;  cq  are  called  eigenvalues  and  <|>,(0)  are  called  eigenfunctions. These 
eigenpairs  ((q,  <j),)  depend  on  the  material  properties,  the  geometry  and  the  boundary  conditions  [<t>;(0) 
are  smooth  vector  functions].  The  A,-  are  coefficients  which  depend  on  the  loading.  Because  of  their 
close  analogy  to  stress  intensity  factors  in  linear  elastic  fracture  mechanics,  At  are  called  generalized 
stress  intensity  factors  (GSIFs).  In  the  case  of  heat  transfer  problems  they  are  called  generalized  flux 
intensity  factors  ( GFIFs). 

The  stresses  in  the  same  neighborhood  can  be  computed  from  the  displacements  given  by  Eq.  (1)  and 
the  material  properties  as: 


{o}  =  £  A{r  1  {^(0)} 

i  =  1 


(2) 
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where  \|t,(0)  depend  on  the  eigenfunctions  <|),(0)  in  Eq.  (1)  and  the  material  coefficients.  It  is  clear 
from  Eq.  (2)  that  when  a,  <  1,  the  stresses  become  singular  for  r=0.  For  additional  details  see  Ref. 
[2]-[5]. 

Three-dimensional  singularities  are  considerably  more  difficult  to  analyze  than  two-dimensional 
ones,  where  only  one  type  of  singularity  exists.  In  3D  in  a  neighborhood  of  the  edges  and  the  vertices 
the  singular  behavior  is  different. 

Edge  Singularities:  If  a  coordinate  system  (x,y,z)  is  located  at  an  edge,  with  the  y-axis  along  the  edge, 
then  there  are  three  edge  GSIFs  which  are  y-dependent:  Aj(y),  Ajrfy)  and  Ajjrfy).  These  edge  GSIFs 
are  analytic  along  the  edge,  however  they  become  themselves  singular  when  this  edge  intersects  with 
a  free  plane,  at  a  vertex.  In  the  neighborhood  of  an  edge- vertex  type  geometry  the  GSIFs  can  be  repre¬ 
sented  once  again  by  vertex  and  vertex-edge  stress  intensity  factors.  For  example,  A,(y)  =  Vs  .  yy,,i  + 
smoother  terms.  j  ’J 

Vertex  Singularities:  In  the  neighborhood  of  a  vertex,  and  away  from  edge-vertex  geometry,  the  dis¬ 
placement  field  can  be  represented  by  only  one  vertex  intensity  factor  and  the  corresponding  eigen- 
pairs.  Investigating  the  mathematical  behavior  of  the  singularities  in  3-D  is  an  active  field  of  research 
in  the  mathematical  community,  and  the  decomposition  of  the  displacement  field  into  singular  and 
regular  terms  is  documented  in  some  recent  papers. 

3.2  Formulation 

The  solution  of  second  order  elliptic  boundary 
value  problems  (BVP)  in  three-dimensions  in  the 
vicinity  of  any  singular  point,  can  be  decomposed 
into  three  different  forms,  depending  whether  the 
singular  point  is  in  the  neighborhood  of  an  edge,  a 
vertex  or  an  intersection  of  the  edge  and  the  ver¬ 
tex.  Mathematical  details  of  the  decomposition  can 
be  found  e.g.  in  Ref.  [8]-[ll]  and  the  references 
therein.  A  representative  three-dimensional 
domain  denoted  by  Q,  which  contains  typical  3D 
singular  points  is  shown  in  Figure  2.  Vertex  singu¬ 
larities  arise  in  the  neighborhood  of  the  vertices  A,-, 
and  the  edge  singularities  arise  in  the  neighbor¬ 
hood  of  the  edges  singularities  A Close  to  the 
vertex/edge  intersection,  vertex-edge  singularities 
arise. 

Consider  the  simplest  elliptic  BVP,  the  Laplace 
equation: 


Ag  A10 


Ar  Vertex  / 

Ajj-  Edge  between  A,  and  Aj 


FIGURE  2.  Typical  3D  singularities. 
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V2«  =  0  in  Q  <3> 

u  =  0  on rD c dQ,  §ii  =  oonrNcaa  (4) 

3  n 

where  u(x j,  xj,  x$)  denotes  the  temperature  field  (in  the  following  x\,  x^  and  x$  will  be  either  Carte¬ 
sian,  cylindrical  or  spherical  coordinates),  and  rD  u  =  3£2 .  It  shall  be  assumed  that  curved  edges 
that  intersect  at  vertices  do  not  exist,  at  that  crack  faces,  if  any,  lie  in  a  plane. 

Edge  Singularities:  We  first  examine  the  edges  denoted  by  Ay  which  connect  the  vertices  A{  and  Ay 
Moving  away  from  the  vertex  a  distance  8/2,  and  creating  a  cylindrical  subdomain  of  radius  r  =  R 

with  the  edge  Ay  as  its  axis,  we  define  a  subdomain  in  the  vicinity  of  the  edge  denoted  bye8  R(Ay) . 

Figure  3  shows  the  edge  singularity  subdomain  e8  R(A,2) .  We  restrict  our  attention  to  domains  having 
straight  edges. 


A, 


FIGURE  3.  The  edge  neighborhood  e8,E(A12). 


The  solution  in  £5  R  can  be  decomposed  as  follows: 


k  s  s 

«(r,0,  z)  =  *s(z)  r“‘( In 0/^(6)  +  v(r,Q,z)  (5) 

*=  If  =  o 
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where  S  >  0  is  an  integer  which  is  zero  unless  ak  is  an  integer,  afc+j  >  ak  are  called  edge  eigenvalues; 
cifofz)  are  called  the  edge  flux  intensity  functions  (EFIFs),  are  analytic  in  z  but  can  become  very  large 
as  they  approach  one  of  the  vertices;  and  fks(Q),  called  eigenfunctions,  are  analytic  in  0.  The  function 
v(r;0,z)  belongs  to  Hg(e),  the  Sobolev  space,  where  q  can  be  as  large  as  required  and  depends  on  K. 
We  shall  assume  that  ak  for  k  <  K  is  not  an  integer,  therefore  Eq.  (5)  becomes: 

K 

M(r,0,  z)  =  ^  ak(z )  r“‘  /*(0)  +  v(r,  0,z)  (6) 

*=  l 

Vertex  Singularities:  A  sphere  of 
radius  p  =  8,  centered  in  the  vertex  A\ 
for  example,  is  constructed  and  inter¬ 
sected  by  the  domain  £1  Then,  a  cone 
having  an  opening  angle  <|>  =  a  is  con¬ 
structed  along  every  edge  intersecting 
at  Aj,  and  removed  from  the  previ¬ 
ously  constructed  subdomain,  as 
shown  in  Figure  4. 

The  resulting  vertex  subdomain  is 
denoted  by  VgfAjl,  and  the  solution  u 
can  be  decomposed  in  Vg (A\)  using  a 
spherical  coordinate  system  by: 


L  Q 

M(P,<(),  0)  =  II  blq  pY'(lnp)\(<t>,0)  +  v(p,<t),0)  (7) 

/  =  i?  =  o 

where  Q  >  0  is  an  integer  which  is  zero  unless  y/  is  an  integer,  y/+1  >  y i  are  called  vertex  eigenvalues, 
and  hiq(i |>,0),  called  the  eigenfunctions,  are  analytic  in  <)>  and  0.  The  b iq  are  called  vertex  flux  intensity 
factors  (VFIFs).  The  function  v(/;0,z)  belongs  to  where  q  depends  on  L.  We  shall  assume  that  Y/ 
for  l  <  L  is  not  an  integer,  therefore,  Eq.  (7)  becomes: 

L 

M(p,<t>,  0)  =  ^ b ,  pY'  h,W,  0)  +  v(p,  <(>,  0)  (8) 

/=  ! 
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Vertex-Edge  Singularities:  The  most  compli¬ 
cated  decomposition  of  the  solution  arises  in  case 
of  vertex-edge  intersections.  For  example,  let  us 
consider  the  neighborhood  where  the  edge  A12 
approaches  the  vertex  Aj.  A  spherical  coordinate 
system  is  located  in  the  vertex  A],  and  a  cone  hav¬ 
ing  an  opening  angle  <J>  =  cr  with  its  vertex  coincid¬ 
ing  with  A  |  is  constructed  with  A12  being  its 
center  axis.  This  cone  is  terminated  by  a  ball¬ 
shaped  basis  having  a  radius  p  =  8,  as  shown  in 
Figure  5. 

The  resulting  vertex-edge  subdomain  is  denoted 
by  V£Se(Aj,  A12),  and  the  solution  u  can  be 
decomposed  in  VeggCAj,  Aj2): 


FIGURE  5.  The  vertex-edge  neighborhood 

^eS,e(Al>  A12)* 


9)  =  jr  jr  ^aksi 


k=ls  =  0\l  =  1 


Y i 

p  +  m 


ks 


(p) 


Q 

1 

l=\q=0 


(sin<t>)“‘l  ln(sin<P)]sg*j(e)  +  II  c/9pY'(lnp)%((t>,e)  +  v(p,<j>,  0)  (9) 


where  mks( p)  is  analytic  in  p;  g^(0)  is  analytic  in  0,  and  t>,0)  is  analytic  in  <|)  and  0.  The  function 
v(r,Q,z)  belongs  to  WfVz.)  where  q  can  be  as  large  as  required  depending  on  L  and  K.  Again  we  shall 
assume  that  y,  for  l  <  L  is  not  an  integer,  and  ak  for  k  <  K  is  not  an  integer,  therefore,  Eq.  (9)  becomes: 


H(p,<)),  0)  = 


*=  1 


^akl  pY'  +  m*(p) 


V- 1 


L, 

(sine)°‘gt(0)  +  y  cf  pY'  /*/(<(>, 0)  +  v(p,((),  0) 


/=  l 


(10) 


The  eigenvalues  and  the  eigenfunctions  are  associated  pairs  (eigenpairs)  which  depend  on  the  mate¬ 
rial  properties,  the  geometry,  and  the  boundary  conditions  in  the  vicinity  of  the  singular  point/edge 
only.  Similarly,  the  solution  for  problems  in  linear  elasticity,  in  the  neighborhood  of  singular  points/ 
edges  is  analogous  to  Eq.  (5)-Eq.  (10),  the  differences  are  that  the  equations  are  in  vector  form  and  the 
eigenpairs  may  be  complex.  For  general  singular  points  the  exact  solution  uEX  is  generally  not  known 
explicitly,  i.e.,  neither  the  exact  eigenpairs  nor  the  exact  EFIFs,  VFIFs  are  known,  therefore  numerical 
approximations  must  be  found. 
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4.1  Task  1 :  Eigenpalrs  for  edge  singularities  in  heat  transfer  problems 

Separation  of  variables  was  utilized  for  the  Laplace  equation  in  the  neighborhood  of  the  edge  sin¬ 
gularity  as  shown  in  Eq.  (6).  Each  eigenpair  r°^4(0)  is  independent  of  z  and  satisfies  the  Laplace 
equation  over  the  plane  (r,0)  which  is  perpendicular  to  the  edge.  This  is  exactly  a  two-dimensional 
problem  and  therefore  the  computation  of  eigenpairs  was  performed  using  the  modified  Steklov 
method  described  in  Ref.  [3]  and  summarized  below. 

The  algorithm  developed  for  two-dimensional  heat-transfer  problems  was  extended  to  three-dimen¬ 
sion  by  considering  a  slice  contained  in  the  plane  normal  to  the  edge  singularity.  The  region  denoted 
by  in  Figure  6,  represents  a  cross  section  along  the  edge  singularity.  The  edge  is  assumed  to  be 


FIGURE  6.  Cross  section  of  an  edge  singularity  for  the  modified  Steklov  formulation. 


along  the  local  z  axis,  and  (r,  0)  are  the  coordinates  of  a  cylindrical  system  located  at  the  intersection 
of  the  normal  plane  and  the  singular  edge.  The  area  around  the  singularity  is  internally  divided  into 
sub-domains. 

This  ‘2D-intemal  mesh’  is  arranged  in  a  circular  ring  around  the  singularity  in  such  a  way  that  the  ele¬ 
ment  boundaries  coincide  with  the  material  interfaces  as  shown  in  Figure  7. 
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Material  1 


Material  3 


Elem.  1 


The  number  of  elements  of  the  internal  mesh 

is  controlled  by  the  number  of  material  inter-  Material  2  Material  1 

faces  around  the  singularity.  The  largest  solid 
angle  for  a  single  element  is  limited  to  120°. 

Therefore,  if  the  partition  of  a  single  material 
is  larger  than  120°,  the  element  is  divided  in  Material  3 
two  or  more  partitions.  For  each  element  of 
the  ‘internal  mesh’  the  corresponding  stiff¬ 
ness  and  trace  matrices  are  computed,  and 
after  assembly,  the  following  eigenvalue 
problem  is  obtained: 

01)  Elem.  3  Sm^ulanty 

where  [K$]  is  the  condensed  stiffness  matrix, 

[M]  is  the  trace  matrix  computed  by  integra-  FIGURE  7-  Typical  ‘internal  mesh’  around  the  singularity. 

tion  on  the  circular  boundary  segments  on 
Q*r,  and  {mr}  is  the  vector  of  coefficients 

that  correspond  to  the  degrees  of  freedom  associated  with  the  circular  boundaries  of  the  internal 
mesh’.  The  solution  of  the  eigen-problem  given  by  the  above  equation  yields  approximation  for  the 
eigenvalues  oq  and  the  corresponding  eigenvectors.  The  steps  and  fundamentals  for  obtaining  the  sys¬ 
tem  described  by  Eq.  (11)  are  described  in  the  following. 

By  formulating  the  weak  form  over  Q*R,  the  singularity  is  excluded  from  the  domain  of  interest  such 
that  the  accuracy  of  the  finite  element  solution  does  not  deteriorate  in  its  vicinity.  On  the  boundaries 
Tj  and  T2  consider  either  zero  temperature  or  zero  flux  boundary  conditions: 


FIGURE  7.  Typical  ‘internal  mesh’  around  the  singularity. 


u  =  0 


9«  _  n 

or  =  0 
an 


on  T/ ,  i  =  1 , 2. 


In  Q*r,  the  temperature  field  u  can  be  represented  as  follows: 


u  =  r7(e) 


Differentiating  Eq.  (13),  and  considering  the  boundary  f3  we  can  write: 


^  =  (a  /R)u 
or 


and  a  similar  condition  on  T^. 

Multiplying  the  Laplace  equation  by  a  test  function  v  in  Hl(Q*R),  integrating  over  the  domain  Q  R 
and  using  Green’s  theorem,  the  modified  Steklov  weak  form  is  obtained: 
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Seek  a €91,  0 *ue  H‘(QR) 

B(u,v)  - (Nr(u,v)  +  Nr.(u,v))  =  a(MR(u,v)  +  MR.(u,v)),  V  ve  h'(Qr)  (15) 

where 

B(u,v)  =  J  J([D]v)r[£]([D]u)</Q,  (16) 

Mr(u,V)  =  \[v[Q}T[E)[Q]u}r  =  R  dB,  (17) 

e 

Nr(u,v)  =  J[v[0]r[£]([D(9)]«)]r  =  *  dB,  (18) 

6 

and  [£>],  [Z/0^],  [Q\  and  [/?]  are  given  as  follows: 


[D]  = 

V 

dx 

[ Q 1  =  cos0  , 

(-sinfljA 

d_ 

[sine 

coseS  _ 

[E]  is  the  2D  material  matrix  reduced  from  the  3D  material  matrix. 

Remark  1.  The  domain  £2*R  does  not  include  singular  edge,  hence  no  special  refinement  of  the  finite 
element  mesh  is  required. 

Remark  2.  The  formulation  of  the  weak  form  was  not  based  on  the  assumption  that  the  material  is 
isotropic,  and  in  fact  can  be  applied  to  multi-material  anisotropic  interface. 

The  domain  Q*R  is  divided  into  finite  elements  through  a  meshing  process,  as  described  before.  The 
polynomial  basis  and  trial  functions,  {'Fj},  are  defined  on  a  standard  element  in  the  T|  space  such 
that  -1<  £  <1,  -1  <  r|  <1.  The  temperature  field  is  then  expanded  in  terms  of  the  known  polynomial 
shape  functions  multiplied  by  a  set  of  unknown  coefficients  {m^}: 

“  =  KJw  <2°) 

The  entries  of  the  unconstrained  stiffness  matrix  corresponding  to  B(u,v)  are  given  by  (see  Ref.  [3]): 

Kij  =  J  J  ([D]{'FI.})r[£][D]{T;.}J£2  (21) 
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Similar  expressions  are  obtained  for  the  matrices  [NR],  [NR*]  and  [MR],  [MR*]  associated  with  the 
bilinear  forms  NR(u,v )  and  MR(u,v). 

Denoting  the  set  of  all  unknown  coefficients  by  {utot},  and  the  set  of  coefficients  associated  with  T3 
and  r4  by  {h^},  the  following  eigen-problem  is  obtained: 


«*]  - 1  Nr]  -  [\.]){«w,}  =  <*([*#,]  +  [^.]){««}  =  a[M]{«^}  (22) 

The  vector  which  represents  the  total  number  of  nodal  values  in  Q*r  can  be  divided  into  two  vectors 
such  that  one  contains  the  coefficients  {uR},  the  other  contains  the  remaining  coefficients:  {utot}T  = 
{  {ur}t,  {w,„}r}.  By  eliminating  {uin},  the  reduced  eigen-problem  is  obtained: 


[Ks]{ur]  =  OL[M]{ur}  (23) 

The  solution  of  the  eigen-problem  given  by  Eq.  (23)  yields  approximations  for  eigenpairs  with  high 
accuracy,  efficiency  and  robustness. 

The  procedure  implemented  in  Stress  Check  for  the  computation  of  the  eigenpairs  requires  the  user 
only  to  click  with  the  mouse  cursor  along  the  singular  edge.  The  program  then  determines  a  cutting 
plane  normal  to  the  singular  edge  at  the  pick  point  and  creates  the  partitions  indicated  in  Figure  7 
based  on  the  material  properties  of  the  elements  intersecting  the  cutting  plane.  Once  these  partitions 
are  identified,  the  elemental  matrices  are  computed  and  the  system  of  equations  is  assembled  and 
solved. 

4.2  Task  2:  Generalized  flux  intensity  factors  for  edge  singularities  in  heat  transfer  problems 

The  algorithm  for  the  computation  of  the  GFIFs  is  based  on  L2  projection  of  the  finite  element  solu¬ 
tion  into  the  space  of  functions  characterized  by  the  asymptotic  expansion  in  terms  of  the  eigenpairs. 
The  algorithm  can  be  summarized  as  follows: 

Consider  the  domain  around  a  singular  edge  for  a  heat  transfer  problem  shown  in  Figure  8.  The  tem¬ 
perature  field  can  be  expanded  around  the  singular  edge  in  terms  of  the  eigenpairs  and  the  generalized 
flux  intensity  factors: 


N 

u(r,  6)  =  ^  A /'/,(&)  =  L4>JM} 

i  =  1 


(24) 


where  At  are  the  GFIFs,  oq  are  the  eigenvalues  and  ^(6)  are  the  corresponding  eigenfunctions,  with 
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Let  Uf^r,Q)  be  the  finite  element  solution  of  the  tem¬ 
perature  distribution  around  the  singular  edge.  Then 
the  L2  projection  of  upp  into  the  domain  QR  is  charac¬ 
terized  by  the  following: 

j  (u-uFE)vdA  =  0,  for  all  ve  S  (25) 

ClR 

where  S  is  the  space  of  eigenfunctions,  and  v(r,0)  is  a 
test  function  given  by: 

N 

v(r,0)  =  ££/'/,(  9)=  iBjW  (26) 

i  =  1 


FIGURE  8.  Typical  cross-section  of  an  edge 
singularity  through  point  R 


Substituting  Eq.  (24)  and  Eq.  (26)  into  Eq.  (25)  and 
rearranging  we  have: 


L*J  f  Wl*]dA{A}  =  |Bj  f  {*}uFEdA  (27) 

Eliminating  B,  from  Eq.  (27)  and  noting  that  dA=rdrdQ,  the  following  system  is  obtained: 

e2  r2  e2 

JJ{<KL<Mr4r40{A}  =  JJ{<t )}uFErdrdQ  (28) 

Rx  0j  /f,  0, 

From  the  definition  of  {  <}>}  in  Eq.  (24)  and  integrating  in  the  radial  direction,  the  system  of  equation 
reduces  to: 


KijAj  =  Rj 


Kij  = 


(a,  +  a,  +  2)  (a,  +  a,  +  2)  e2 

*■  2,a"l 


I  J 


R,  02 


Rj  =  j  jr(a‘+Ufj(Q)uFEdrde 


r,  e, 


(29) 


Solving  the  system  of  equations  represented  by  Eq.  (29)  over  the  domain  gives  the  values  of  A,- 
(the  GFIFs).  The  L2  projection  provides  excellent  approximation  when  the  generalized  intensity  fac- 
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tor  along  the  singular  edge  is  either  constant  or  linear,  but  the  approximation  deteriorates  for  cases  in 
which  the  variation  is  of  higher  order.  Nevertheless  the  current  implementation  provides  a  very  robust 
algorithm  for  a  large  class  of  practical  problems. 

The  procedure  implemented  in  Stress  Check  for  the  computation  of  the  generalized  flux  intensity  fac¬ 
tors  requires  the  user  to  first  solve  the  heat  transfer  problem  given  the  topology,  material  properties 
and  boundary  conditions.  After  the  finite  element  solution  is  available,  the  post-processing  operation 
requires  only  to  click  with  the  mouse  cursor  along  the  singular  edge.  The  program  then  determines  a 
cutting  plane  normal  to  the  singular  edge  at  the  pick  point  and  computes  the  eigenvalues  and  corre¬ 
sponding  eigenvectors  as  explained  in  the  previous  section.  With  the  eigenpairs  and  the  temperature 
distribution  obtained  from  the  finite  element  solution  the  program  constructs  the  system  given  by  Eq. 
(29),  the  solution  of  which  provides  the  GFIFs  as  shown  in  the  next  section. 

4.3  Task  3:  Model  problems  -  Eigenpairs  and  generalized  flux  intensity  factors 

The  procedures  described  in  the  previous  sections  for  the  computation  of  the  eigenpairs  and  the 
corresponding  generalized  flux  intensity  factors  for  edge  singularities  in  steady  state  heat  transfer 
(scalar)  problems,  are  illustrated  in  the  following. 

Scalar  problem  1:  Clamped-free  crack  in  isotropic  material.  Circular  domain  of  unit  radius  with  a 
crack  along  the  positive  x-axis.  One  face  of  the  crack  (Tj)  has  zero  temperature  boundary  condition 
and  the  other  face  (r2)  is  flux  free.  The  outside  boundary  (TR)  has  an  imposed  flux  (Figure  9).  The 
finite  element  mesh  consisting  on  16  solid  elements  is  also  shown  in  Figure  9. 
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The  elements  are  graded  towards  the  singular  edge  in  geometric  projection  with  a  common  factor  of 
0.15. 

The  exact  solution  for  this  problem  is  given  by  (Ref.  [1]): 


«(r,0)  =  -1.35812r1/4sin(0/4)  +  0.970087  r3/4  sin(30/4)  +  ... 


(30) 


The  exact  values  of  the  first  two  eigenvalues  for  this  problem  are:  aj=l/4  and  a2=3/4,  and  the  corre¬ 
sponding  generalized  flux  intensity  factors  are:  A]=l. 358097,  A2=0.970087. 

The  finite  element  solution  was  obtained  for  polynomial  orders  ranging  from  1  to  8.  The  estimated 
relative  error  in  energy  norm  for  the  sequence  of  finite  element  solutions  is  shown  in  Figure  10. 


(1)  Solution  =  SOL,  runs  #1  to  #8 


Run  # 

DOF 

Total  Potential  Energy 

Convergence  Rate 

%  Error 

1 

32 

-2 . 032896794902e+000 

0.00 

31.85 

2 

112 

-2 -228519967775e+000 

0.76 

12.24 

3 

200 

-2 .248067253504e+000 

0.74 

7.96 

4 

352 

-2  -252250724991e+000 

0.31 

6.70 

5 

572 

-2.254400378353e+000 

0.25 

5.95 

6 

876 

-2  .255778943319e+000 

0.22 

5.41 

7 

1280 

-2.256771999437e+000 

0.21 

4.99 

8 

1800 

-2  .257536919443e+000 

0.21 

4.64 

Estimated  Limit  -2 .262400478095123e+000 


FIGURE  10.  Estimated  relative  error  in  energy  norm  for  scalar  problem  1 . 


The  first  and  second  eigenvalues  computed  in  Stress  Check  and  the  corresponding  values  of  the 
GFIFs  are  shown  in  Figure  1 1 .  The  tabular  data  shows  the  convergence  of  the  Aj  and  A2  as  a  function 
of  the  number  of  degrees  of  freedom  (DOF)  corresponding  to  eight  finite  element  solutions  obtained 
for  polynomial  orders  ranging  from  1  to  8.As  can  be  seen,  the  numerically  computed  eigenvalues  are 
practically  identical  to  the  theoretical  values.  The  generalized  flux  intensity  factor  are  also  very  close 
to  the  theoretical  values.  The  estimated  limits  shown  in  the  table  are  computed  by  a  projection  to  infi¬ 
nite  number  of  degrees  of  freedom  of  the  three  values  with  the  highest  DOF.  The  value  in  brackets 
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represents  the  percent  deviation  between  the  estimated  limit  and  the  value  corresponding  to  run  #8. 


Number  of  E-pairs:  2,  Solid  angle:  3 .60000Qe+002 
Global  coord,  of  point  along  edge:  X=0.Q,  Y=0.0,  2=0.0 


No.  1,  Eval  =  2.500000e-001 
No.  2,  Eval  =  7.500000e-001 


GFIFs 

Isotropic  clamped-free  crack 

(0)  Solution  =  SOL,  runs  1  to  8  (nodes=16-16,angle=360 .0) , 
Generalized  Flux  Intensity  Factors,  Int.  Radius  =0.5 


Run 

DOF 

Radius 

A1 

A2 

1 

32 

5 .00000e-001 

-1.12165e+000 

8 . 40991e-001 

2 

112 

5 .00000e-0Ql 

-1.28522e+000 

9 . 60927e-001 

3 

200 

5 .OOOOOe-OOl 

-1.32626e+000 

9 . 70949e-001 

4 

352 

5 .00000e-001 

-1.33414e+000 

9.70400e-001 

5 

572 

5. OOOOOe-OOl 

-1.33718e+0Q0 

9.70066e-001 

6 

876 

5. OOOOOe-OOl 

-1.34012e+0Q0 

9.70035e-0Ql 

7 

1280 

5. OOOOOe-OOl 

-1.34257e+000 

9.70081e-001 

8 

1800 

5. OOOOOe-OOl 

-1.34437e+000 

9 . 70092e-001 

Est. Limit  Al=-1.353007e+000  (  0.64%) 
Est .Limit  A2=9 .700747e-001  (  0.00%) 


FIGURE  11.  Eigenvalues  and  Generalized  Flux  Intensity  Factors  for  scalar  problem  1. 


Finally,  Figure  12  shows  the  numerically  computed  eigenfunctions  associated  with  the  first  and  sec¬ 
ond  eigenvalues. 

E-function  associated  with  E-value:  2.50000Qe-001  E-function  associated  with  E^/alue:  7.5Q00Q0e-001 


FIGURE  12.  Eigenfunctions  associated  with  the  first  and  second  eigenvalues  for  scalar  problem  1. 
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Scalar  problem  2:  Anisotropic  reentrant  corner.  This  is  a  heat  transfer  problem  in  an  anisotropic 
material  whose  boundary  consists  of  a  90°  reentrant  comer  generated  by  two  flux-free  faces  T j  and 
r2,  which  meet  along  an  edge  as  shown  in  Figure  13.  The  cylindrical  boundary  of  the  domain  (TR)  is 
loaded  by  flux  boundary  condition  which  corresponds  to  the  first  symmetric  eigenfunction  of  the 
asymptotic  expansion  of  u(x,y,z)  about  the  reentrant  edge  as  given  in  Ref.  [6].  The  finite  element  mesh 
created  in  Stress  Check  consisting  of  three  solid  elements  is  also  shown  in  the  figure. 


FIGURE  13.  Anisotropic  reentrant 
mesh. 


Material  properties: 

he c=4.0 
kyy=kzz=  1 .0 


Boundary  Conditions: 

h(0,0,0)  =  0.0 
du/dQ  =  Oon  Tj,  r2 


corner.  Notation  and 


The  exact  values  of  the  first  and  second  eigenvalues  and  the  corresponding  GFIFs  are  aj=2/3,  Aj=1.0 
and  a2=4/3,  A2=0.0,  respectively.  The  finite  element  solution  was  obtained  for  polynomial  orders 
ranging  from  1  to  8.  The  estimated  relative  error  in  energy  norm  for  the  sequence  of  solutions  is 
shown  in  Figure  14.  The  estimated  error  is  around  3%  for  the  solution  with  433  degrees  of  freedom 
(DOF). 


Error  Estimate 

Anisotropic  reentrant  corner 
(0)  Solution  -  SOL,  runs  #1  to  #8 


Run  # 

DOF 

Total  Potential  Energy 

Convergence  Rate 

%  Error 

1 

9 

-8 . 823  659296528e-001 

0.00 

29.06 

2 

28 

-9 . 42592296 1277e-001 

0.59 

14.81 

3 

53 

-9 . 540787 15 1703 e-OOl 

0.61 

10.01 

4 

91 

-9 . 5848 17509 602 e-OOl 

0.56 

7.38 

5 

145 

-9 . 605587799983 e-001 

0.54 

5.74 

6 

218 

-9. 61669155 6281e-001 

0.53 

4.63 

7 

313 

-9 . 623 188417919e-001 

0.52 

3.83 

8 

433 

-9 . 627258673 192e-001 

0.52 

3.24 

Estimated  Limit  -9 . 637358320546628e-001 


FIGURE  14.  Estimated  realtive  error  in  energy  norm  for  scalar  problem  2. 
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The  first  and  second  eigenvalues  computed  in  Stress  Check  and  the  corresponding  eigenvectors  are 
shown  in  Figure  15,  while  the  values  of  the  GFIFs  are  shown  in  Figure  16.  The  tabular  data  shows  the 


E-function  associated  with  E-value:  6.666667e-001  E-function  associated  with  E-value:  1 ,333333e+000 


FIGURE  15.  Eigenfunctions  associated  with  the  first  and  second  eigenvalues  for  scalar  problem  2. 


convergence  of  the  A\  and  A2  as  a  function  of  the  number  of  degrees  of  freedom  (DOF)  correspond¬ 
ing  to  eight  finite  element  solutions  obtained  for  polynomial  orders  ranging  from  1  to  8.  As  can  be 

Mnber  of  E-pairs:  2,  Solid  angle:  2 .700000e+002 
Global  coord,  of  point  along  edge:  3M. 0,  Y-0.0,  Z— 0.147912 
No.  1,  Eval  -  6 . 666667e-001 
No.  2 ,  Eral  -  1 .333333e+000 


Fracture  Extraction 
Anisotropic  reentrant  corner 
(0)  Solution  -  SOL,  runs  1  to  8  <nodes-l-l , angle-270 . 0) , 
Generalized  Flux  Intensity  Factors,  Int.  Radius  -  0.9 


Run 

DOF 

Radius 

Al 

A2 

1 

9 

9 . 000e-~001 

-8 . 850e-001 

-5 . 038e-002 

2 

28 

9 . 000e-001 

-9 . 879e-001 

-2.4430-003 

3 

53 

9 . 000e-001 

-9 . 837C-001 

-4.930e-004 

4 

9X 

9 . 000e-001 

-9 . 951e-001 

6 . 555e-004 

5 

145 

9 . 000e-001 

-9 . 962e-001 

2 . llle-004 

6 

218 

9 . 000e-001 

-9.970e-001 

4 . 884e-005 

7 

313 

9. 000e-001 

-9 . 987e-001 

2.482e-004 

8 

433 

9 . 000e-001 

-9 . 984e-001 

7 • 792e-006 

Est. Limit  Al— 9 . 979279e-001 
Est. Limit  A2-1 . 103098e-004 


FIGURE  16.  Eigenvalues  and  Generalized  Flux  Intensity  Factors  for  scalar  problem  2. 


18  of  68 


Stress/Failure  Analysis  Software  for  Multi-Material  Interfaces 


Implementation  and  results 


seen,  the  numerically  computed  eigenvalues  are  practically  identical  to  the  theoretical  values.  The 
generalized  flux  intensity  factor  are  also  very  close  to  the  theoretical  values. 

Scalar  problem  3:  Two  material  interface.  Two  materials  perfectly  bonded  along  a  common  edge 
satisfying  the  following  equation: 


£,V2m  =  0  in  Q(  (31) 

with  the  following  flux  conditions  along  the  external  boundary  (Figure  17): 

—.  =  k^r1'  '  h^Q)  +  a.^*2  '/i2(e)]  on  I>  3Q,  i  =  1,2  (32) 

The  material  coefficients  are:  fcj=10,  the  eigenvalues  are:  aj=0.73 1691779,  a2=  1.268308221, 
and: 


/*,(©)  = 


*2(0)  = 


{ 

{ 


cos [( 1  —  #)0]  +  Cj  sin [( 1  —  o)0]  — >  on  Tj 
Cj cos [(1  -a)0]  +  c2c3sin[(l  -a)0]  on  T2 

cos[(l  +  <3)0]-c3sin[(l  +  a)0]  on  Tj 
CjCOs[(l  +  a)0]-c2c3sin[(l  +  a)0]  — >  on  T2 


(33) 


(34) 


where  c1=6. 31818181818182,  c2=-2. 68181818181818,  c3=0. 64757612580273,  and 
a=0.26830822 130025.  The  exact  solution  for  this  problem  is  (Ref.  [7]): 

u(r,6)  =  Alra'h{(Q)  +  A2ra2h2(Q)  (35) 


where  Aj  =  A2=  1. 

A  sequence  of  finite  element  solutions  was  obtained  for  the  mesh  shown  in  Figure  17  consisting  of  6 
solid  elements  graded  in  geometric  progression  towards  the  singular  edge  with  a  common  factor  of 
0.15.  The  estimated  relative  error  in  energy  norm  for  this  problem  is  shown  in  Figure  18.  As  can  be 
seen  from  the  results,  the  global  error  of  the  solution  is  well  under  1%  for  p=8  (716  DOF). 
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Material  properties: 

On  Qj.-  kxx  =kyy=kzz=l0.0 
On  Gig:  kxx  =kyy=kzz=  1.0 


Boundary  Conditions: 

du/d r  =f(r,Q)  on  rl5  T2 


FIGURE  17.  Two  material  interface.  Notation  and  mesh. 


Error  Estimate 
Two  material  Interface 
(0)  Solution  =  SOL,  runs  #1  to  #8 


Run  # 

DOF 

Total  Potential  Energy 

Convergence  Rate 

%  Error 

1 

13 

-5.135011591076e+001 

0.00 

73.67 

2 

44 

-1.110656827166e+002 

1.60 

10.51 

3 

81 

-1.115029819039e+002 

0.36 

8.46 

4 

142 

-1.122893439218e+002 

3.39 

1.26 

5 

230 

- 1. 123 0 12918 676e+0 02 

1.14 

0.73 

6 

351 

- 1. 123 0 448372 05e+0 02 

0.90 

0.50 

7 

511 

- 1. 123 05482250 0e+0 02 

0.59 

0.40 

8 

716 

-1.123060704408e+002 

0.59 

0.33 

Estimated  Limit  -1.123072750529331e+002 


FIGURE  18.  Estimated  relative  error  in  energy  norm  for  scalar  problem  3. 
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The  first  and  second  eigenpairs  are  shown  in  Figure  19.  Note  that  the  maximum  of  these  functions  are 
located  at  225  degrees  for  the  first  eigenfunction  and  at  154  degrees  for  the  second  eigenfunction.  In 


E-function  associated  with  E-value:  7.31 691 8e-001 


both  cases  the  maximum  is  unity.  For  the  exact  eigenvectors  given  in  Eq.  (33)  for  the  first  eigenvalue 
and  in  Eq.  (34)  for  the  second  eigenvalue,  the  maximum  occur  at  the  same  angular  positions  (225  and 
154  degrees  respectively)  but  their  magnitude  is  different  (6.5525  for  A,j  and  h^). 

The  first  and  second  eigenvalues  and  the  corresponding  flux  intensity  factors  computed  numerically 
are  shown  in  Figure  20.  As  before,  the  approximation  of  the  eigenvalues  is  very  good,  and  the  values 
of  A  j  and  A2  are  also  close  when  considering  the  difference  in  the  amplitude  of  the  eigenvectors  men¬ 
tioned  above.  Referring  to  Eq.  (35),  it  is  apparent  that  the  approximate  and  exact  solutions  give  the 
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same  temperature  field  in  the  neighborhood  of  the  singular  edge.  That  is  because  the  product  of  the 
GSEFs  and  the  eigenvectors  is  the  same. 


Number  of  B-pairs:  2,  Solid  angle:  3 .600000e+002 

Global  coord,  of  point  along  edge:  X=0.0,  Y=0.0,  Z=-5.04279e-002 

No.  1,  Eval  =  7.316918e-001 

No.  2,  Eval  =  1.268308e+000 


Fracture  Extraction 
Two  material  Interface 

(0)  Solution  =  SOL,  runs  1  to  8  (nodes=17-17,angle=360.0) , 
Generalized  Flux  Intensity  Factors,  Int.  Radius  =0.5 


Run 

DOF 

Radius 

A1 

A2 

1 

13 

5 .000e-001 

5.379e+000 

8.777e-001 

2 

44 

5 .000e-001 

6.40Se+000 

6.875e+000 

3 

81 

5 .000e-001 

6.536e+000 

6.883e+000 

4 

142 

5 .000e-001 

6.545e+000 

6.586e+000 

5 

230 

5 .000e-001 

6.544e+000 

6.549e+000 

6 

351 

5 .000e-001 

6.546e+000 

6.551e+000 

7 

511 

5 .000e-001 

6.547e+000 

6 .552e+000 

8 

716 

5 .000e-001 

6.548e+000 

6.552e+000 

Est. Limit  Al=6.549912e+000  (  0.03%) 
Est. Limit  A2=6.552459e+000  (  0.00%) 


FIGURE  20.  Eigenvalues  and  Generalized  Flux  Intensity  Factors  for  scalar  problem  3. 
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Finally,  Figure  21  shows  the  convergence  graph  of  and  A2  as  a  function  of  the  number  of  degrees 


Two  materia!  Interface 

(0)  Solution  -  SOL  runs  1  to  8  (nodes»1 7-1 7,angle-360.0), 
Generalized  Flux  Intensity  Factors,  Ini  Radius  *  0.5 


Legend 

— * —  A1 
- A2 


FIGURE  21.  Convergence  of  A1  and  A2  for  scalar  problem  3. 


of  freedom  (DOF).  It  is  apparent  from  the  results,  that  the  generalized  flux  intensity  factors  are  practi¬ 
cally  independent  of  the  discretization. 


Scalar  problem  4:  Isotropic  reentrant  corner.  To  illustrate  the  quality  of  approximation  when  the 
generalized  flux  intensity  factor  varies  along  the  singular  edge,  consider  the  reentrant  comer  shown  in 
Figure  22. 
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This  is  a  heat  transfer  problem  of  an  isotropic  material  whose  boundary  consists  of  a  90°  reentrant 
comer  generated  by  two  flux-free  faces  Tj  and  T2,  which  meet  along  an  edge.  The  cylindrical  bound¬ 
ary  of  the  domain  (TR)  and  the  front  (z=L)  and  back  (z=0)  faces  have  Diritchlet  boundary  conditions, 
that  is  the  temperature  distribution  w(r,0,z)  is  imposed  using  the  expansion: 

u(r,Q,z)  =  cosia^)  +  A2(z)ra2cos(a2Q) 


A,  (z)  =  an+a12z,  A2(z)  =  a21+a22z 


where  a,=/7t/G)  are  the  exact  eigenvalues  and  the  variation  of  the  GFIFs  A\  and  A2  along  the  singular 
edge  is  linear.  Selecting  a j  j  =  a2j  =  1.0  and  aj2  =  a22  =  0.5,  the  exact  values  of  the  GFIFs  along  the 
singular  edge  are: 


A,(z)=  A2(z)=  Aex  =  1.0  +  0.5  z 

A  sequence  of  finite  element  solutions  was  obtained  for  the  mesh  shown  in  Figure  22  consisting  of  9 
solid  elements  graded  in  geometric  progression  towards  the  singular  edge  with  a  common  factor  of 
0.15.  The  estimated  relative  error  in  energy  norm  is  shown  in  Figure  23.  With  this  mesh,  the  error  at 


Error  Estimate 
Isotropic  reentrant  corner 
(0)  Solution  ■  SOL,  runs  #1  to  U 8 


Run  U 

DOF 

Total  Potential  Energy 

Convergence  Rate 

H  Error 

2 

9 

1 . 102985683990e+001 

0.00 

11.12 

3 

16 

1. 10740952  2  03  Oe-fOOl 

0.28 

9.16 

4 

45 

1. 108501491406e+001 

0.07 

8.61 

5 

93 

1. 1109444852 18e+001 

0.24 

7.23 

6 

171 

1.1153823  672 17e+0Q 1 

1.17 

3.54 

7 

288 

1. 1164866123 63 e+001 

1.49 

1.63 

8 

453 

1. 116705996903 e+001 

1.49 

0.83 

Estimated  Limit  1. 116782792307470e+0Ql 


FIGURE  23.  Estimated  realtive  error  in  energy  norm  for  scalar  problem  4. 


p=8  (453  DOF)  is  less  than  1%.  The  corresponding  GFIF  values  obtained  using  the  l ?  projection 
algorithm  implemented  in  Stress  Check  are  shown  in  Table  1  at  various  positions  along  the  singular 
edge.  As  can  be  seen,  the  numerical  results  are  practically  identical  to  the  exact  values. 


24  of  68 


Stress/Failure  Analysis  Software  for  Multi-Material  Interfaces 


Implementation  and  results 


TABLE  1.  Numerical  and  exact  values  of  GFIFs  along  singular  edge 


z 

Ai(z) 

A2(z) 

aex 

0.053 

1.026 

1.026 

1.0265 

0.664 

1.332 

1.332 

1.3320 

1.333 

1.667 

1.667 

1.6665 

1.967 

1.983 

1.983 

1.9835 

4.4  Task  4:  Eigenpairs  for  edge  singularities  in  elasticity  problems 

The  elastostatic  displacement  field  in  three-dimensions  in  the  vicinity  of  an  edge  singularity  can 
be  decomposed  in  terms  of  the  eigenpairs  and  the  corresponding  stress  intensity  functions.  Mathemat¬ 
ical  details  on  the  decomposition  are  found  in  Ref.  [8],  [10]  and  the  references  therein.  Elastic  edge 
singularities  associated  with  anisotropic  materials  and  multi-material  interfaces  have  been  less  inves¬ 
tigated,  however.  Analytical  methods  as  in  Ref.  [12],  [13]  provide  procedures  for  the  computation  of 
eigenpairs,  but  require  very  involved  mathematical  operations. 

In  the  following,  a  description  of  the  numerical  procedure  implemented  in  Stress  Check  for  the  com¬ 
putation  of  eigenpairs  for  edge  singularities  in  elasticity  is  addressed.  In  the  neighborhood  of  a  typical 
edge  singularity  the  displacement  fields  can  be  written  as: 

K  s 

u(r,Q,  z)  =  ££atj(z)  rat(lnr)X(Q)  +  w(r,Q,z)  (36) 

*=  li  =  0 

where  S  >  0  is  an  integer  which  is  zero  for  most  practical  problems,  and  aks(z)  are  analytic  in  z  and  are 
called  the  generalized  stress  intensity  functions.  The  vector  function  w(r,Q,z )  belongs  to  [/^(eg  R)]3  , 
where  q  depends  on  K.  Considering  the  case  S  =  0,  Eq.  (36)  becomes: 

K 

u(r,Q,z)  =  ^ak(z)  fk(6)  +  w(r, 0,z)  (37) 

k  =  i 

The  tractions  on  the  boundaries  are  denoted  by  T  =  (Tx  Ty  tz)t,  and  the  Cartesian  stress  tensor  by 
a  =  ax  ov  oz  %xy  iyi  tzx  .  Body  forces  are  not  considered  in  the  vicinity  of  the  singular  edge. 
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FIGURE  24.  The  modified  Steklov  domain  £lR. 


For  computing  the  eigenpairs,  a  two-dimensional 
sub-domain  Q,R  is  constructed  in  a  plane  perpen¬ 
dicular  to  the  edge  bounded  by  r=R  and  r=R*  as 
shown  in  Figure  24.  On  the  boundary  0=0  and 
0=£D12  °f  the  sub-domain  QR  homogenous  trac¬ 
tions,  homogeneous  displacement  constraints  or  a 
combinations  of  these  are  prescribed. 

In  view  of  Eq.  (37)  the  displacement  vector  in  QR 
has  the  following  functional  representation: 


u  =  a(z)M/v(6)  =  a(z)rXf(0) 


and  the  in-plane  variation  of  the  displacements  is 
written  as: 


u(r,  0)  =  -j— 

a(z ) 


Following  the  steps  presented  in  detail  in  Ref. 
[14],  an  eigenvalue  problem  is  cast  in  a  weak  form  over  a  two  dimensional  domain  involving  all  three 
displacement  components: 

SeekaeC,  0*ue  [//' (£2^) ] suchthat,  Vve  [//'(Q*)]3 


B(u,  v)-[Nr(u,  v)-Nr.(u,  v)  1  =  X[Mr(u,v)-Mr,(u,  v)] 


T 

B(u,v)  =  f  j“,2|([Ar)ar+[/te]^)v|  [£]|([Ar]ar+[Ae]^)«|r4r4e 


CO  1 2 

^fi(M,v)=f  (v)r[Af]r[£][Ae]39«|  £ 

~  ~  J  n  ~  ~\r  =  R 


Mr(u,  v)  =  J  '\v)TlAr]T[E][Ar]u\r  =  RdQ 
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where  [E]  is  a  6x6  material  matrix  with  up  to  21  independent  coefficient  for  a  completely  anisotropic 
material;  d^=  d/dr,  3q=  3/30;  and 


COS0 

0 

0 

-sin0 

0 

0 

0 

sin0 

0 

0 

COS0 

0 

0 

0 

0 

[Ae]  = 

0 

0 

0 

sin0  cos0 

0 

COS0 

-sin0 

0 

0 

0 

sin0 

0 

0 

COS0 

0 

0 

COS0 

0 

0 

-sin0 

Remark  1.  The  test  and  trial  functions  have  three  components,  but  the  domain  over  which  the  weak 
eigen-formulation  is  defined  is  two-dimensional  and  excludes  any  singularity.  Therefore  the  applica¬ 
tion  of  the  p- version  of  the  finite  element  method  (FEM)  for  solving  Eq.  (41)  is  very  efficient. 

Remark  2.  The  bilinear  forms  NR  and  NR*  are  non-symmetric  with  respect  to  «  and  v  thus  they  are 
not  self-adjoint.  As  a  consequence  the  ‘minimum  principle’  does  not  hold,  and  any  approximation  of 
the  eigenvalues  obtained  using  a  finite  dimension  subspace  of  [//‘(Q*)]3  cannot  be  considered  as  an 
upper  bound  of  the  exact  ones,  and  the  monotonic  behavior  of  the  error  is  lost  as  well.  Convergence  is 
assured  under  a  general  proof  provided  in  Ref.  [16],  however. 

Next  consider  the  discretization  of  the  weak 
formulation  given  by  Eq.  (40)  using  the  p- 
version  FEM  over  a  finite  dimensional  sub¬ 
space  of  [//'(Q*)]3.  Assuming  that  the 
domain  ClR  consist  of  three  different  materi¬ 
als  as  shown  in  Figure  25,  let  us  divide  it 
into  three  finite  elements  as  shown.  Let  us 
consider  a  typical  element  bounded  by  0j 
and  02  (element  1  in  Figure  25).  A  standard 
element  is  considered  in  the  £,  T)  plane  such 
that  -1<  £  <1,  -1  <  T|  <1  over  which  the  poly¬ 
nomial  basis  and  trial  functions  are  defined. 

These  standard  elements  are  mapped  into  the 
real  elements  by  proper  mapping  functions 
(for  details,  see  Ref.  [15]  chapters  5  and  6). 

The  functions  £  and  v  are  expressed  in 
terms  of  the  basis  functions  4/j-(£,T|)  in  the 
standard  plane: 

V,  ...  Vw  0  ...  0  0  ...  0  c, 

u  =  o  ...  0  V,  ...  Vat  0  ...  0  ...  =  1VF]{C}  (43) 

0  ...  0  0  ...  0  V,  ...  VaJLc3n_ 
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where  C,  are  the  amplitudes  of  the  basis  functions  and  \|/,  are  the  products  of  the  integral  of  the  Leg¬ 
endre  polynomials  in  £,  and  T|.  Since  u  and  v  lie  in  the  same  space,  we  define  similarly  v  =  [T]{B}. 

The  unconstrained  stiffness  matrix  corresponding  to  the  bilinear  form  B  in  Eq.  (41)  is  given  by: 

[K]  =  J62j(jArR+  [Ae]y)m  J  [E]j([ArR  +  [Ae]^)m  jr4r40  (44) 

Next  we  define  the  matrices: 


tfP]  = 


-PjSinG  ... 

-P^sinG  0 

0 

0 

0 

0 

0  P,cosG  . 

..  PN  COS0 

0 

0 

0 

0  0 

0 

0 

0 

P,cos0  ... 

/  / 

PN cos0  -P,  sin0  . 

..  -PwsinG 

0 

0 

0 

0  0 

0 

P,COS0  ... 

PN COS0 

0 

0  0 

0 

-P,sinG  ... 

-P^sinG 

[PI  = 


P,COS0  . 

.  P^cosG 

0 

0 

0 

0 

0 

0 

PjSinG  ... 

P^sinG 

0 

0 

0  . 

0 

0  ... 

0 

0  ... 

0 

/>1sin0  . 

.  PN sinG  P{c os0  ... 

P^cosG 

0 

0 

0 

.  0 

0  ... 

0 

PjSinG  ... 

P^sinG 

0  . 

0 

0  ... 

0 

PjCosG  ... 

PN  cosG 

where  P(<£)  for  i>3  are  the  integrals  of  the  Legendre  polynomials,  and  P1(^)=(l-^)/2,  P2(0=(1+^V2- 
Therefore,  the  bilinear  form  N  can  be  written  as: 


Nr(u,v)  =  [B]t  |  [pf[£][3P]  |„=_,4  1 C ]  =  [Bfl^llC] 


The  entries  in  [NR]  are  computed  using  Gauss  quadrature.  Similarly,  the  expression  for  the  bilinear 
form  M  is  evaluated  as  follows: 

MrCu,v)  =  IPlfeV  [PflE][~P]  L.^llCl  =  [Bf[MR][C]  (46) 

-i  J 

The  matrices  [NR*]  and  [M^*]  have  the  same  values  as  those  of  [NR]  and  [MR],  but  of  opposite  sign. 
This  is  because  the  shape  functions  on  boundaries  T3  and  T4  are  the  same  (except  for  a  sign  change), 
and  so  is  the  mapping  to  the  standard  element.  Denoting  the  set  of  amplitudes  of  the  basis  functions 
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associated  with  the  boundary  T3  by  [Cy,  and  those  associated  with  the  boundary  r4  by  [C#*],  the 
eigenpairs  can  be  obtained  by  solving  the  following  generalized  eigenvalue  problem: 

[*][C]  -  [Nr][Cr]  +  [AW[C*.]  =  XCEA/^HC^]  -  [M*,]  [C*,])  (47) 

Augmenting  the  coefficients  of  the  basis  functions  associated  with  T3  with  those  associated  with  r4 
and  denoting  them  by  the  vector  [C/^*],  Eq.  (47)  becomes: 

[*][C]-[A^.][Crr.]  =  X[Mrr.][Crr.]  (48) 

The  vector  that  represents  the  total  number  of  nodal  variables  in  Q.R  is  divided  into  two  vectors,  one 
containing  the  variables  [C^*]  and  the  other  the  rest  of  the  variables:  [C]T  =  {[C^*]r,  [C,n]r}.  By 
partitioning  the  stiffness  matrix  [A(],  Eq.  (48)  can  be  written  as: 

\k]-[Nrr.]  [Krr._J  _  ^["[^^*1  toil  *j 

[*J  J[cinJ  [[0]  [0]J[cinJ 

The  relation  in  Eq.  (49)  is  used  to  eliminate  [Cjn]  by  static  condensation,  thus  obtaining  the  reduced 
eigenvalue  problem: 

[KjlfCflyj*)  =  (50) 

where 

[*si  = 


In  Eq.  (50)  [Ay  is  a  full  matrix.  However  the  order  of  the  matrices  is  relatively  small,  so  the  solution 
of  the  eigenvalue  problem  (using  Cholesky  factorization  to  compute  [Afy1)  is  not  expensive.  Note 
that  the  derivation  was  performed  for  a  single  element.  For  multiple  elements  along  T3  and  T4  the  for¬ 
mulation  is  identical  and  the  matrices  [AT],  [NR]  and  [MR\  are  obtained  by  an  assembly  procedure. 

4.5  Task  5:  Generalized  stress  intensity  factors  for  edge  singularities  in  elasticity  problems 

As  in  the  case  of  the  generalized  flux  intensity  factors,  the  algorithm  for  the  computation  of  the  GSIFs 
is  based  on  a  L2  projection  of  the  finite  element  solution  into  the  space  of  functions  characterized  by 
the  asymptotic  expansion  in  terms  of  the  eigenpairs.  The  algorithm  can  be  summarized  as  follows: 
Consider  a  section  perpendicular  to  a  singular  edge  of  an  elasticity  problem  as  shown  in  Figure  26. 
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The  displacement  field  corresponding  to  the  exact  solution 
can  be  expanded  around  the  singular  edge  at  point  P  in  terms 
of  the  eigenpairs  and  the  generalized  stress  intensity  factors: 

ux(r,  0)  AT 

u(r,Q)  =  -  uy(r,  0)  >  =  ^4,/'/(0)  (51) 

.  «f(r,0)J  '  =  1 

where  At  are  the  GSIFs,  X,-  are  the  eigenvalues  and 
{/(0)}={/x(0X/j(0)>/z(9)}  are  the  corresponding  eigenfunc¬ 
tions.  Let  {uF£^r,Q)}  be  the  finite  element  solution  of  the  dis¬ 
placement  field  around  the  singular  edge.  Then  the  L2 
projection  of  { up^}  into  the  domain  Qr  is  characterized  by  the 

FIGURE  26.  Typical  cross-section  of  an  following: 
edge  singularity  through 
point  P. 

J  (m  -  uFE)vdA  =  0,  for  all  v  €  S  (52) 

Qr 

where  S  is  the  space  of  eigenfunctions,  and  {v(r,0)J,  is  given  by: 

N 

v(rfi)  =  £fi/£(0) 

I  =  1 

Substituting  {«}  and  {v}  into  Eq.  (52)  and  rearranging  we  have: 

L«J  J  ({<M|AJ  +  {^iLfvJ  +  {<t>z}L<t>J)<M{/i}  =  LS_|J({<UkF£|;c  +  {<t>,  }*F£,,  + (53) 

where: 

(<U  =  rX‘{fx},  {<|>v}  =  rK{fy},  {♦,}  =  rX‘{fz}  (54) 

Eliminating  Bt  from  Eq.  (53)  and  noting  that  dA=rdrdQ,  the  following  system  is  obtained: 

*2  e2  *2  e2 

e,  Rt  e, 

From  the  definition  of  {((>}  in  Eq.  (54)  and  integrating  in  the  radial  direction,  the  system  of  equation 
reduces  to: 
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¥;  =  RJ 


(Xj  +  Xj  +  2)  (X(  +  X;  +  2)  0. 


Ku  =  ~  2  +  x!‘x - J  ^(0)4(0)  +/v,(0)/v/Q)  +/z/(6)/yfe)]^e 

'  7  ®i 

r2  e2 

^/  =  J  Jr  (fijUFE\x  +fyjUFE\y  +fzjUFE\z^r^ 


Solving  the  system  of  equations  represented  by  Eq.  (56)  over  the  domain  QR  gives  the  values  of  A, 
(the  GSIFs).  Note  that  only  a  line  integral  is  required  to  compute  Ky  and  an  area  integral  to  compute 
Rj.  The  displacement  vector  from  the  finite  element  solution  {wpg}  has  to  be  evaluated  relative  to  the 
local  coordinate  system  located  at  the  extraction  point  P. 

The  procedure  implemented  in  Stress  Check  for  the  computation  of  the  generalized  stress  intensity 
factors  is  the  same  as  described  for  the  heat  transfer  problem.  First  the  elasticity  problem  is  solved  for 
the  given  topology,  material  properties  and  boundary  conditions.  After  the  finite  element  solution  is 
available,  the  post-processing  operation  requires  only  to  click  with  the  mouse  cursor  along  the  singu¬ 
lar  edge.  The  program  then  determines  a  cutting  plane  normal  to  the  singular  edge  at  the  pick  point 
and  computes  the  eigenvalues  and  corresponding  eigenvectors  as  explained  in  the  previous  section. 
With  the  eigenpairs  and  the  displacement  field  obtained  from  the  finite  element  solution  the  program 
constructs  the  system  given  by  Eq.  (56),  the  solution  of  which  provides  the  GSIFs  as  illustrated  in  the 
next  section. 

Stress  Intensity  Factors.  Utilizing  most  of  the  functions  and  procedures  developed  for  the  computa¬ 
tion  of  GFIFs/GSIFs,  an  algorithm  was  implemented  to  compute  the  stress  intensity  factors  for  cracks 
in  isotropic  materials  using  the  contour  integral  method  (CIM).  The  mode  1  and  2  stress  intensity  fac¬ 
tor  are  computed  from  the  finite  element  solution  at  any  position  along  the  crack  front  using  the  CIM 
as  described  in  Ref.  [15],  page  227.  The  procedure  is  outlined  in  the  following. 

For  plane  stress  and  plane  strain,  the  stress  intensity  factors  in  mode  I  and  mode  II  are  given  by: 


Kj  =  JlK/ 


=  *Jln/ 


where  A-^m\  m  =  1,  2,  is  the  first  term  in  the  asymptotic  expansion  of  the  solution  in  the  neighbor¬ 
hood  of  the  crack  tip.  Let  (jt,y)  represent  a  local  Cartesian  coordinate  system  located  at  any  point 
along  the  crack  front  in  such  a  way  that  the  local  z-axis  is  tangent  to  the  crack  front  at  that  point,  and 
the  x-axis  is  in  the  direction  of  the  crack  face.  Let  T  be  a  circle  of  radius  p  centered  at  the  point  along 
the  crack  front,  and  assume  that  p  is  sufficiently  close  to  the  crack  (Figure  27).  It  has  been  shown  that 
Aj  can  be  computed  from  a  line  integral  as: 
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,("0 


i 


w  r„„  —  «__r 

m  FE  FE 


SWj)ds 


(58) 


where  Wm  is  an  extraction  function,  TFE  is  the  trac¬ 
tion  vector  along  T  computed  from  the  finite  ele¬ 
ment  solution  uFE,  and  T^Wm)  is  the  traction  vector 
along  r  due  to  the  extraction  function.  The  traction 
and  displacement  vectors  must  be  computed  in  rela¬ 
tion  to  the  local  coordinate  system  attached  to  the 
crack  front. 


The  terms  in  Eq.  (58)  are  computed  from  the  finite  FIGURE  27.  Contour  around  crack  front, 

element  solution  of  the  cracked  body  along  the  con¬ 
tour  T  as  follows: 


TFE  - 


TJx,  y) 


l_ Ty(x,  y) 


JFE 


a  cos0  +  t  sine 

A  AJ' 

LVCOS0  +  OySinejf£ 


UFE  ~ 


ux(x,  y ) 


uy(x,  y) 


(59) 


JFE 


The  terms  in  Eq.  (58)  due  to  the  extraction  functions  are  given  by: 


■=  {'F(m>(0)})  T(Wm)  =  -^-{T<m)(0)} 


(60) 


D 


where  G  is  the  shear  modulus,  and  D(,)  =  7t( 2k  -  1),  Z)(2)  =  7C(2k  +  3),  where  k  depends  on  the  Pois¬ 
son’s  ratio  (v).  For  plane  strain  k  =  3  -  4v,  and  for  plane  stress  k  =  (3  -  v)/(l  +  v).  Also,  the  vector 
functions  that  depend  on  the  variable  0  are  given  by: 


(  net  39 

Ik  —  cos  -  -  -  cos— 

l  2)  2  2  2 

(  A  .  8  1  .  30 

lK+2jSm2  2SI"  2 


{XP(2)(0» 


(  3^  .  0  1  .  30 

k  +  -  sm-  +  -sin— 

l  2j  2  2  2 

(  y\  0  1  30 


(61) 
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(3  0  1  50^  . 

[2COS2  +  2COStJCOS0  + 

fjsin^-jsin^lcosfi  +  l 
\2  2  2  2) 


'\  .  56  i  .  ev  a 

-sin—  --sin-  sin0 
,2  2  2  2) 

5  0  1  50A  .  . 

2cos2_2cosy)sine 


(  7  ■  0  1  •  50^  Q 
—  sin- --sin—  cos0  + 

l  2  2  2  2) 

(l  .  50  3  0)  .  / 

-sin  —  +  -  cos-  cos0  + 

\2  2  2  2/  l 


1  50  3  0)  .  0 

-cos—  +  -COS-  SH10 

2  2  2  2) 

1.0  1  .  50V 

2  2  2  2  ) 


(62) 


In  general,  a  three-dimensional  problem  is  neither  a  plane  stress  nor  a  plane  strain  problem.  Therefore 
the  question  is  how  to  determine  the  value  of  K  for  the  extraction  function  in  Eq.  (61).  The  procedure 
to  determine  which  value  of  K  to  use  is  based  on  computing  the  stress  components  ox,  oy,  az  at  a  point 
in  front  of  the  crack  (see  Figure  27).  Next  define  y  as  the  ratio: 

Y  =  —2s. . 

0,  +  °, 


If  y  is  close  to  Poisson’s  ratio  (v),  then  the  condition  approaches  that  of  plane  strain.  If  is  close  to  zero, 
then  is  plane  stress.  The  approach  implemented  in  Stress  Check  was  to  compute  the  value  of  y  at  a 
point  along  the  integration  path  (r  =  p,  0  =  0),  and  apply  the  following  rule  to  determine  k: 

3-4v  if  (y>  0.85v) 
if  (y<  0.85v) 

1  +v 


The  advantages  of  this  implementation  are  improved  computational  speed  and  the  ability  to  include 
curved  crack  fronts.  The  improvements  in  performance  are  realized  because  the  eigenvalues  and  cor¬ 
responding  eigenvectors  need  not  be  computed  prior  to  computing  the  stress  intensity  factors.  An 
example  of  the  implementation  is  presented  in  the  next  section. 

4.6  Task  6:  Model  problems  -  Eigenpairs  and  generalized  stress  intensity  factors 

The  procedures  described  in  the  previous  sections  for  the  computation  of  the  eigenpairs  and  the 
corresponding  generalized  stress  intensity  factors  for  edge  singularities  for  elasticity  problems  in  3D, 
are  illustrated  in  the  following. 

Elasticity  problem  1:  Isotropic  reentrant  corner.  The  first  example  of  the  implementation  of  the 
algorithm  for  the  computation  of  the  eigenpairs  and  generalized  stress  intensity  factors  for  edge  sin- 
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gularities,  consider  the  traction-free  isotropic  L-shaped  domain  shown  in  Figure  28  loaded  along  the 
boundaries  by  the  Mode  1  and  Mode  2  stress  components  as  given  in  Ref.  [15].  The  load  is  applied  to 
the  faces  labeled  1,  through  4  in  Figure  28.  The  two  free  faces  perpendicular  to  the  singular  edge  were 
constrained  from  motion  in  the  z-direction. 

The  finite  element  mesh  created  in  Stress  Check  consisting  of  18  solid  elements  with  two  layers  of 
geometrically  graded  elements  towards  the  reentrant  comer  is  also  shown  in  the  figure.  The  elements 
were  graded  with  a  common  factor  of  0.15.  The  problem  was  solved  for  polynomial  order  ranging 


from  p=l  to  8.  The  estimated  relative  error  in  energy  norm  for  the  sequence  of  eight  solutions  is 
shown  in  Figure  29.  For  p=8  (5338  DOF)  the  global  error  of  the  finite  element  solution  is  less  than 


1%. 


ErrorEstimate 
L-SHAPED  DOMAIN 

(0)  Solution  =  SOL,  runs  #1  to  #8 


Run  # 

DOF 

Total  Potential  Energy 

Convergence  Rate 

%  Error 

1 

85 

-4.299854922149e-001 

0.00 

24.66 

2 

307 

-4.54B143062914e-001 

0.87 

8.11 

3 

553 

-4.571877574981e-001 

1.32 

3.73 

4 

988 

-4.57 6467 179110e- 001 

1.10 

1.97 

5 

1630 

-4.577446575364e-001 

0.80 

1.32 

6 

2533 

-4.5777 80333023e-001 

0.61 

1.01 

7 

3751 

-4.5 779 41449 00 8e-001 

0.54 

0.82 

6 

5338 

-4.5 7 8 038073 468e- 001 

0.54 

0.67 

Estimated  Limit  -4.5782461Q9968524e-G01 


FIGURE  29.  Estimated  relative  error  in  energy  norm  for  elasticity  problem  1. 
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The  exact  values  of  the  first  three  eigenvalues  and  the  corresponding  GFIFs  are  A,j=0.544483, 
Ai=1.0,  A^=2/3,  A2=0.0  and  A,3=0.908529,  A3=1.0,  respectively.  The  eigenvalues  computed  in  Stress 
Check  are  shown  in  Figure  30,  together  with  the  three  components  of  the  first  and  third  eigenvectors. 
Note  that  all  the  eigenvalues  are  real. 


REAL  E-functions  associated  with  E-vaiue:  5.444837e-001  ♦  1m  0.000000e+000 


1.0 

0.8 

0.6 

0.4 

0.2 

0.0 

-0.2 

-0.4 

-0.6 

-0.8 


100  200  300 

Theta 


FIGURE  30.  Eigenpairs  for  elasticity 
problem  1. 


Number  of  E-pairs:  4,  Solid  angle:  2.700000e-K302 
Global  coord,  of  point  along  edge: 

X=0 .00000e+000,  Y=0 .00000e+000,  Z=-4.94063e-005 
No.  1,  Eval=  5.444837e-001  +  I  0 . DOOOOOe+QOO 

No.  2,  Eval=  6.666667e-001  +  I  0 .000000e+000 

No.  3/  Eval=  9 .085292e-001  +  I  0 .OOOOOOe+OOO 

No.  4,  Eval=  1.000000e+000  +  I  0 .OOOOOOe+OOO 

Legend 
*  lb<J 
UyJ 

— + —  UzJ 


REAL  EHiinctions  assodaled  with  E-value:  9.085292e-001  +  Im  0.000000e+000 


2.0 


4: 


Legend 

-*■■■  Ux  3 

*  (JyJ 

— * —  Uz_3 


-2.0 


The  values  of  the  GSIFs  computed  by  the  numerical  procedure  described  in  the  previous  section  are 
shown  in  the  Figure  31.  The  tabular  data  shows  the  GSIFs  computed  from  each  finite  element  solu¬ 
tion.  The  graph  shows  the  values  of  A1  and  A3  as  a  function  of  the  number  of  degrees  of  freedom 
(DOF).  The  stress  intensity  factors  are  practically  coincident  with  the  exact  values. 
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GSIFS 

L- SHAPED  DOMAIN 

(0)  Solution  =  SOL,  runs  1  to  8  (nodes=22-22,angle=270 .0) , 
Generalized  Stress  Intensity  Factors,  Int.  Radius  =0.5 


Run 

DOF 

Radius 

A1 

A2 

A3 

A4 

1 

85 

5 -OOOOOOe-OOl 

-7.107936e-001 

1 . 446139e-G14 

-9.391302e-001 

8.031775e-001 

2 

307 

5 .OOOOOOe-OOl 

-9 .38203 9e-001 

1.445316e-006 

-9 -890786e-001 

8.060001e-001 

3 

553 

5.000000e-001 

-9.836161e-001 

3 ,560110e-008 

-9.983398e-001 

8.104599e-001 

4 

988 

5.000000e-001 

-9.926223e-001 

3 .474569e-011 

-1.00Q214e+000 

8 ,112030e-001 

5 

1630 

5 .000000e-0Ql 

-9.936874e-001 

4.783306e-010 

-1.00Q123e+000 

8.114920e-001 

6 

2533 

5 .OOOOQOe-DOl 

-9 .946648e-001 

4.768257e-010 

-9 .999739e-001 

8.113539e-001 

7 

3751 

5 .000000e-001 

-9 .955081e-001 

-1.794964e-008 

-9.999516e-001 

8 ,113397e-001 

6 

5338 

5 .000000e-001 

-9.961335e-001 

“9 .7549G8e-Q10 

-9.999505e-001 

0 .1133Q3e-001 
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Elasticity  problem  2:  Crack  in  isotropic  material.  Consider  the  case  of  a  through-thickness  crack  in 
a  thick  isotropic  plate  under  tension  loading  as  shown  in  Figure  32.  After  the  finite  element  solution  is 


FIGURE  32.  Through-thickness  crack  in  thick  plate.  Undeformed  and  deformed  configurations. 

obtained,  the  procedure  for  the  computation  of  the  stress  intensity  factors  requires  the  user  simply  to 
point  and  click  with  the  mouse  cursor  to  the  desired  location  along  the  crack  front. 

The  procedure  implemented  in  Stress  Check  determines  a  cutting  plane  normal  to  the  tangent  to  the 
crack  edge  at  the  pick  point  and  extract  the  global  components  of  the  stresses  and  displacements  along 
a  circular  path  contained  in  the  cutting  plane.  These  stresses  and  displacements  are  projected  into  the 
cutting  plane  and  integrated  with  the  extraction  function  to  compute  the  contour  integral  described  in 
the  previous  section. 

Figure  33  shows  the  results  of  the  stress  intensity  factors  (SEFs)  obtained  from  a  set  of  finite  element 
solutions  for  polynomial  orders  ranging  from  p=l  to  6  at  three  locations  along  the  crack  front.  Z=0.5 
represents  the  center  of  the  plate  and  Z=0.0  is  one  of  the  free  faces.  Note  the  variation  of  the  SIFs 
along  the  crack  front,  with  the  largest  value  of  K1  at  the  center  of  the  plate.  The  estimated  limits  for 
the  SIFs  are  also  included  which  indicate  the  strong  convergence  characteristics  of  the  implemented 
extraction  procedure. 

Elasticity  problem  3:  Inclusion  problem.  Consider  a  composite  body  consisting  of  two  dissimilar 
isotropic,  homogeneous  and  elastic  wedges,  perfectly  bonded  along  the  interfaces.  The  body  is  loaded 
on  the  circular  boundary  by  the  stress  field  corresponding  to  the  exact  solution  of  the  asymptotic 
expansion  about  the  singular  point  as  given  in  Ref.  [17].  The  normal  displacement  is  set  to  zero  in  the 
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Contour  Integral  Method:  JLngle-3  60,  Radius-2 .50e-01 

Edge 

location:  X-8. 

OSe-17,  Y«1.21e-17 

/ 

Z-S.00e-01  1 

Run  # 

DOF 

K1 

K l  - 

1 

X08 

2 . 344084e+000 

8 . 085629e-002 

2 

3  69 

4 . 301251e+000 

3 

.570169e-002 

3 

678 

4. 861055e+000 

3 

. 683432e-002 

4 

lies 

5. 160177e+000 

3 ,276876e-002 

5 

1914 

5.268396e+000 

3 

.567037e-002 

6 

2913 

|S.300811e+000  | 

3 

.498816e-002 

1  Estimated  Limits 

Limit 

value  for  K1 

...  5. 319954e+000 

(diff:  0.36%) 

Limit 

value  for  K2 

...  3 . 497732e-002 

(diff:  0.03%) 

Edge 

location:  X«8.05e-17,  Y-1.21e-17 

9 

Z*2 . Sle-01  | 

Run  U 

DOF 

K1 

K2 

1 

108 

2.352391e+000 

*\ 

\861335e-002 

2 

3  69 

4 . 200244e+000 

3 

. 580445e-002 

3 

678 

4 . 802555e+Q00 

3 

.711076e-002 

4 

1185 

5.153232e+000 

3 

.338650e-002 

5 

1914 

5.198848e+000 

3 

. 617087e-002 

6 

2913 

| 5. 189015e+000  | 

3 

.574476e-002 

1  Estimated  Limits 

Limit 

value  for  K1 

...  5 . 146324e+000 

(diff:  0.83%) 

Limit 

value  for  K2 

...  3 . 5S8739e-002 

(diff:  0.44%) 

Edge 

location:  X-8.0Se-17,  Y-1.21e-17 

/ 

Z»2.64e-02  | 

Run  # 

DOF 

K1 

K2 

1 

108 

2. 359888e+000 

*7 

. 6S8902e-002 

2 

3  69 

3 .935169e+000 

3 

.6091236-002 

3 

678 

4 . 649050e+000 

3 

.783665e-002 

4 

1185 

4. 899167e+000 

3 

.462604e-002 

5 

1914 

4.893078e+000 

3 

.787856e-002 

6 

2913 

|4.862767e+000  | 

3 

. 679010e-002 

1  Estimated  Limits 

Limit 

value  for  K1 

...  4.856227e+000 

(diff:  0.13%) 

Limit 

value  for  K2 

...  3 . 68242 Se-002 

(diff;  0.09%) 

FIGURE  33.  Stress  intensity  factors  for 
elasticity  problem  2. 


two  lateral  faces.  Figure  34  shows  the  description  of  the  problem  and  the  finite  element  mesh  consist¬ 
ing  of  10  solid  elements.  The  first  three  eigenvalues  characterizing  the  stress  singularity  along  the 
edge  are:  X,  =  0.5124722,  ^  =  0.6001175  and  X3  =  0.7309757. 

A  sequence  of  finite  element  solutions  was  obtained  for  polynomial  orders  ranging  from  1  to  8.  The 
estimated  relative  error  in  energy  norm  is  shown  in  Figure  35.  The  first  four  eigenvalues  computed  in 
Stress  Check  are  shown  in  Figure  36,  together  with  the  corresponding  generalized  stress  intensity  fac¬ 
tors  (A,-,  /  =  1  to  4).  Note  that  all  the  eigenvalues  are  practically  identical  to  the  analytical  values,  and 
the  excellent  convergence  characteristics  of  the  GSIFs  with  increasing  number  of  degrees  of  freedom. 

Finally,  the  three  components  of  the  first  three  eigenvectors  are  shown  in  Figure  37. 
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Error 

Inclusion  problem 
(1)  Solution  ■  SOL,  runs  #1  to  #8 


Run  M 

dof 

Total  Potential  Energy 

Convergence  Rate 

*  Error 

1 

41 

-5 . 453568540785e-004 

0.00 

31.96 

2 

154 

-5 . 980512 13  6567e-004 

0.72 

12.40 

3 

287 

-6 . 041965203489e-004 

0.86 

7.24 

4 

520 

-6 . 0577598 15888e-004 

0.58 

5.14 

5 

668 

-6 . 063591919594e-004 

0.44 

4.11 

6  . 

1361 

-6. 0666902593 lSe-004 

0.40 

3.43 

7 

2029 

-6. 06860912 5881e-004 

0.39 

2.93 

8 

2902 

-6 . 06988 602 82 95e-004 

0.39 

2.55 

Estimated  Limit  -6 . 073828017459567e-004 


FIGURE  35.  Estimated  relative  error  in  energy  norm  for  elasticity  problem  3. 


Elasticity  problem  4:  Angles-ply  anisotropic  laminate.  Consider  the  singularities  associated  with  a 
composite  laminate  with  ply  properties  typical  of  a  high-modulus  graphite  epoxy  material.  The  mate¬ 
rial  coefficients  in  the  principal  direction  of  the  fibers  are: 

i?£=20xl06  psi,  £'7=2.1x106  psi,  Gif=Gjf=0 . 8 5 x 1 06  psi,  V£j=V77=0.21 
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Number  of  E-pairs:  4,  Solid  angle:  3.600e+002 
Global  coord,  along  edge:  X-0.0,  Y-0.0,  Z-7.45e-03 
No.  1,  Eval*  5. 124722e-001  +  I  Q.OOOOOQe+OOO 

No.  2,  Eval-  6 . 001175e-001  +  X  0.000000e+000 

No.  3/  Eval-  7.309757e-001  +  I  0.000000e+000 

No.  4,  Eval-  1.000000e+000  +  I  0 . 000000e+000 


GSir 

Inclusion  problem 

(2)  Solution  -  SOL,  runs  1  to  8  (nodes-7-7, angle-3 60. 0) , 
Generalised  Stress  Intensity  Factors,  Int.  Radius  -  0.2 


Run 

DOF 

Radius 

A1 

JL2 

A3 

A4 

1 

41 

2 .000e-001 

8. 566e-005 

-1 . 135e-016 

3 . 971e-004 

2 . 663e-004 

2 

154 

2 .000e-001 

1.286e-004 

-1.684e-007 

3 ,952e-004 

1 . 704e-004 

3 

287 

2 .000e-001 

1.423e-004 

-8 . 186e-009 

4.023e-004 

1 . 694e-004 

4 

520 

2. OOOe-OOl 

1.466e-004 

-8. 658e-012 

4.046e-004 

1.724e-004 

5 

868 

2. OOOe-OOl 

1.490e-004 

-4 . 625e-Q12 

4.055e-004 

1.728e-004 

6 

1361 

2. OOOe-OOl 

1.505e-004 

-4. 364e-012 

4.059e-0Q4 

1.733e-004 

7 

2029 

2. OOOe-OOl 

1.515e-004 

4 . 630e-012 

4.062e-004 

1 . 73  6e-004 

8 

2902 

2. OOOe-OOl 

1.522e-004 

4.145e-012 

4.063e-004 

1.738e-004 

Est. Limit  Al-1.551560e-004  (  1.91%) 


FIGURE  36.  Eigenvalues  and  GSIFs  for  elasticity  problem  3. 


where  L  and  T  refer  to  the  fiber  and  transverse  directions,  respectively.  There  are  three  plies  in  the 
laminate  and  they  are  oriented  at  +  or  —  45  degrees  about  the  global  y-axis  (Figure  38).  Of  interest  are 
the  eigenvalues  associated  with  the  free  edge  singularity  between  two  plies,  and  those  associated  with 
a  delamination  (edge  crack)  between  two  plies.  In  the  case  of  the  free-edge  singularities  the  eigenval¬ 
ues  reported  in  Ref.  [18]  are: 

X1  =  0.974424,  7^2,3  =  1-88147  ±  i  0.234005,  X.4>5  =  2.511526  ±  i  0.792817 
For  the  delamination  (edge  crack),  the  reported  values  are: 

X12  =  0.50  ±  i  0.0343,  A,3  =  0.50 

For  the  computation  of  eigenvalues  in  Stress  Check  it  is  necessary  to  perform  at  least  one  solution.  In 
this  case  the  mesh  shown  in  Figure  38  consisting  of  6  solid  elements  was  solved  for  p=l  only.  The 
eigenvalues  were  then  extracted  at  the  free  edge  singularity  and  at  the  crack  singularities  indicated  in 
Figure  38.  The  results  are  shown  in  Figure  39.  As  can  be  seen  they  match  perfectly  the  values 
reported  in  Ref.  [18]. 

Elasticity  problem  5:  Composite  patch.  Consider  the  a  composite  patch  bonded  to  a  metallic  struc¬ 
ture  as  shown  in  Figure  40.  Assuming  that  the  thickness  of  the  adhesive  is  negligible,  the  composite 
patch  is  in  full  contact  with  the  metal.  Therefore  we  want  to  investigate  the  edge  singularity  along  the 
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REAL  E-funclions  associated  with  E-value:  5.124722e-001  +  Im  0.000000e+000 


REAL  E-functions  associated  with  E-value:  6.0011 75e-001  +  Im  0.000000e+000 


-100  100  200  300 

Theta 


Legend 


FIGURE  37.  Eigenvectors  for  elasticity 
problem  3. 


composite  patch  as  a  function  of  the  termination  angle  a.  The  material  coefficients  in  the  principal 
direction  of  the  fibers  are: 

El=2 OxlO6  psi,  2i7=2.1xl06  psi,  Gij<=Gff=0.S5xl06  psi,  \Lj=Vjj=0.21 

with  the  fibers  rotated  an  angle  (3  about  the  y-axis.  For  the  aluminum  plate  the  following  properties 
are  considered:  is^lO.SxlO6  psi  and  v=0.3.  We  will  investigate  the  eigenpairs  obtained  with  a  ply  ori¬ 
entation  of  (3=45°  while  changing  the  termination  angle  from  a  from  20°  to  90°.  The  finite  element 
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Singular  edge 


-45 

+45 

-45 


Delamination  (edge  crack) 

FIGURE  38.  Elasticity  problem  4.  Notation  and  Mesh. 


Number 

of  E-pairs:  6,  Solid  angle:  160 

Global 

X- 

coord,  of  point  along  edge: 

1.0,  Y-0.2,  Z-0.S18353 

No. 

1, 

Eval- 

9.744236e-001 

+ 

I  0 . 000000e+000 

No. 

2, 

Eval- 

1 . 000000e+000 

+ 

I  0 . 000000e+000 

No. 

3, 

Eval- 

l.OOOOQOe+OOO 

+ 

I  0 . 000000e+000 

No. 

4, 

Eval- 

1.8B1466e+000 

+ 

I  2 . 340048e-001 

No. 

5, 

Eval- 

1 .861466e+000 

+ 

I-2.340048e-001 

No. 

6, 

Eval- 

1 . 999999e+000 

+ 

I  0 . 000000e+000 

No. 

7, 

Eval- 

2.511439e+000 

+ 

I  7 . 929259e-001 

No. 

0, 

Eval- 

2.511439e+000 

+ 

I-7.929259e-001 

Edge  Singularity 


Delamination 


Number  of  E-pairs:  4,  Solid  angle:  360 
Global  coord,  of  point  along  edge: 
X-0.5,  Y-0.1,  Z-0. 468659 


No.  1,  Eval- 
No.  2,  Eval- 
No .  3 ,  Eval- 
No.  4/  Eval- 


5. 000003e-001  +  I  3 . 434452e-002 
5.000003e-001  +  1-3 ,434452e-002 
5 . OOOOlle-ODl  +  I  0 . 000000e+000 
1 . 000000e+000  +  I  0 . 000000e+000 


FIGURE  39.  Eigenvalues  for  elasticity  problem  4. 
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FIGURE  40.  Elasticity  problem  5.  Composite  patch  over  aluminum  plate. 

mesh  consisting  of  three  solid  elements  used  for  the  numerical  analysis  is  shown  in  Figure  41.  Since 


FIGURE  41.  Finite  element  mesh  fo  elasticity  problem  5. 


the  objective  of  the  analysis  is  the  computation  of  eigenpairs,  it  is  sufficient  to  perform  a  finite  ele¬ 
ment  analysis  at  p=l.  The  results  of  the  analysis  using  the  procedures  implemented  in  Stress  Check 
are  shown  in  Figure  42.  Note  that  the  first  and  second  eigenvalues  become  smaller  as  the  termination 
angle  increases  and  their  magnitudes  are  minimum  (Xj  =0.664 18,  A,2=0.750854)  for  a  =  90°,  which 
corresponds  to  a  solid  angle  of  270°.  If  the  patch  would  have  been  made  of  aluminum  and  the  termi¬ 
nation  angle  a  =  90°,  the  first  two  eigenvalues  would  be:  X]  =0.544837,  ^2=0. 666667,  which  is  a 
‘stronger’  edge  singularity.  These  results  are  in  agreement  with  those  presented  in  Ref.  [18].  Finally, 
Figure  43  shows  the  three  components  of  the  eigenvector  associated  with  X]  =0.664 18  for  a  =  90°. 
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Number  of  E-pairs 
No. 


No. 


Eval- 

Eval- 


2,  So /id  angle: 
+  I 


8.283 671e-001 
9. 661601e-001 


+  I 


200 

0. 000000e+000 
0 . 000000e+000 


Number  of  E-pairs:  2,  Solid  angle:  210 
No.  1,  Eval-  7.866517e-001  +  I  0.000000e+000 
No.  2,  Eval-  9 . 466995e-001  +  I  0.000000e+000 


Number  of  E-pairs:  2,  Solid  angle:  220 
No.  1,  Eval-  7.5B6947e-001  +  I  0 . 000000e+000 
No.  2,  Eval-  9.234772e-001  +  I  0 . 000000e+000 

Number  of  E-pairs:  2,  Solid  angle:  230 
No.  1, 

No.  2, 


Number  of  E-pairs:  2,  Solid  angle:  240 
No.  1,  Eval-  7. 195816e-001  +  I  0 . 000000e+000 

No.  2,  Eval-  8 • 597589e-001  +  X  O.OQOOOOe+OOO 

Number  of  E-pairs:  2,  Solid  angle:  250 
No.  1,  Eval-  7 . 023354e-001  +  I  0.000000e+000 

No.  2,  Eval-  8.206522e-001  +  I  0.000QQ0e+000 


Eval-  7.375793e-001 
Eval-  8.947598e-001 


+  I  0 . OOOOOOe+OOO 


+  I  0. OOOOOOe+OOO 


Number  of  E-pairs:  2,  Solid  angle:  260 
No.  1,  Eval-  6.641136e-001  +  I  0 . OOOOOOe+OOO 

No.  2,  Eval-  7.827444e-001  +  I  0. OOOOOOe+OOO 

Number  of  E-pairs:  2,  Solid  angle:  270 
No.  1,  Eval-  6 . 641812e- 

No.  2,  Eval-  7. 50854  le- 


-001  +  X  0. OOOOOOe+OOO 
■001  +  I  0. OOOOOOe+OOO 


FIGURE  42.  Eigenvalues  for  elasticity  problem  5. 


4.7  Task  7:  Average  stress/strain 

A  set  of  failure  criteria  for  composite-laminated  bonded  structures  based  on  stress  averaging  and 
material  nonlinear  analysis  were  identified  as  a  necessary  complement  to  the  set  of  tools  available  for 
analysis.  These  criteria  include: 

•  Criteria  for  failure  in  the  composite  adherents  (laminate):  Two  failure  modes,  in-plane  failure  and 
out-of-plane  failure,  are  considered  separately  for  the  composite  adherents.  For  in-plane  failure, 
maximum  stress,  maximum  strain,  Tsai-Wu,  Tasi-Hill  and  Hashin  criteria  are  widely  used  in  the 
industry.  For  the  out-of-plane  failure,  a  quadratic  criterion  for  interlaminar  normal  and  shear 
stresses  is  considered.  Because  of  the  stress  singularity  and  high  stress  gradient  near  the  edge  of  the 
laminate  or  near  the  edge  of  a  bonded  joint,  the  use  of  the  average  stress  is  necessary  as  follows: 
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REAL  E-functions  associated  with  E-value:  6.641 81 2e-001  *  Im  0.000000e+000 


Legend 


UxJ 
Uy  1 
UzJ 


FIGURE  43.  Eigenfunctions  associated  with  ^=0.66418  for  a  =  90°.  Elasticity  problem  5. 


where 


(63) 


(64) 


are  the  average  interlaminar  normal  and  shear  stresses  in  the  composite,  a§  is  a  characteristic 
length  (of  the  order  of  the  lamina  thickness),  Gq  is  the  flatwise  tensile  strength  of  the  composite  and 
t0  is  the  interlaminar  shear  strength  of  the  composite.  A  similar  relation  can  be  obtained  for  three 
dimensions,  where  the  line  integral  is  replaced  by  an  area  integral. 

•  Criterion  for  failure  at  the  adhesive/adherent  interface'.  Failure  along  the  interface  is  similar  to  the 
interlaminar  or  out-of-plane  failure  of  the  laminate. 

•  Criterion  for  failure  in  the  adhesive  layer:  Two  major  approaches  have  been  applied  to  determine 
failure  in  the  adhesive  layer:  Stress/strain  based  criteria  and  fracture  or  energy  based  criteria.  Based 
on  the  data  available  and  based  on  the  anticipated  users  of  the  analysis  tool,  a  stress/strain  based 
criterion  is  recommended.  Similar  to  the  interlaminar  stresses,  very  steep  stress/strain  gradients  in 
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the  adhesive  layer  can  be  expected  at  the  end  of  the  bonded  joint.  To  relieve  the  effects  of  this 
stress  concentration,  an  average  strain  criterion  is  therefore  recommended.  This  criterion  can  be 
expressed  as: 


where  the  shear  and  peel  strain  terms  are  averaged  over  a  characteristic  length  a0  of  the  order  of  the 
adhesive  thickness,  and  are  computed  from  a  material  nonlinear  analysis  of  the  bonded  joint.  The 
adhesive  failure  strains,  Eq  and  y0  include  the  plasticity  component  and  are  highly  dependent  on  the 
material  model  used  in  the  analysis. 

A  capability  was  implemented  in  Stress  Check  to  compute  average  stress/strain  along  element  edges, 
element  faces,  element  volumes  and  arbitrary  curves.  The  average  is  understood  in  the  integral  sense, 
and  when  appropriate,  the  user  controls  the  number  of  points  used  for  the  numerical  integration.  If  the 
average  is  computed  along  an  element  edge,  face  or  volume,  then  the  integration  is  performed  using 
Gauss  quadrature,  and  the  number  of  integration  points  is  selected  internally  by  the  program.  If  the 
average  is  performed  along  an  arbitrary  path  which  runs  inside  one  or  more  elements,  then  the  inte¬ 
gration  is  performed  using  the  trapezoidal  rule  and  the  user  controls  the  number  of  integration  points. 

For  example,  if  the  average  stress  az  is  requested  along  an  element  edge  of  length  a0,  the  program 
performs  the  numerical  integration  indicated  in  Eq.  (64).  The  arc  length  a0  is  determined  automati¬ 
cally  based  on  the  screen  selection.  If  the  selected  object  is  an  element  face,  the  integral  is  performed 
over  the  area  of  the  face.  If  the  selected  object  is  an  element,  then  a  volume  integral  is  performed.  The 
length,  area  or  volume  of  the  selected  object  is  also  reported  with  the  corresponding  average  quantity. 
The  procedure  to  compute  the  average  is  illustrated  with  an  example. 

Example:  Consider  a  typical  composite  bonded  joint  under  axial  loading  shown  in  Figure  44.  After 
the  linear  solution  is  obtained,  the  user  selects  Results  >  Points  and  completes  the  extraction  form  for 
the  desired  average  quantity  as  shown  in  Figure  45.  In  this  example,  the  average  of  ay  (labeled  Sy  in 
the  form)  is  computed  along  the  indicated  element  edge  for  a  sequence  of  finite  element  solutions 
ranging  from  p-level=l  to  8.  The  result  shown  includes  the  average  for  each  run  with  the  convergence 
information  and  the  length  of  the  selected  edge. 

The  quantities  for  which  the  average  can  be  computed  include  the  displacement  components  in 
the  global  system,  all  directional  strain  and  stress  components,  principal  strains  and  stresses.  Addi¬ 
tionally,  the  average  can  be  computed  for  quantities  extracted  in  a  local  coordinate  system  and  any 
other  fimction  that  can  be  described  using  the  formula  option  in  Stress  Check.  This  allows  formulat¬ 
ing  an  expression  of  “margin  of  safety”  based  on  average  quantities  for  any  failure  mode  in  the  adhe¬ 
sive,  adherent  or  adhesive/adherent  interface.  For  example,  an  expression  for  the  “margin  of  safety” 
(. MS)  for  the  interlaminar  out-of-plane  failure  of  composite  adherents  can  be  written  as: 
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MS  = 


(66) 


In  this  case,  the  average  of  the  normal  and  shear  stresses  are  computed  first  before  the  expression  of 
MS  is  evaluated. 

4.8  Task  8:  Limit  load  analysis 

The  evaluation  of  a  composite  bonded  joint  requires  a  material  nonlinear  analysis  for  the  design 
load  followed  by  the  determination  of  the  failure  loads  based  on  some  predefined  criteria.  The  evalu¬ 
ation  of  failure  load  by  an  automatic  procedure,  that  can  be  utilized  from  the  handbook  framework  of 
Stress  Check,  required  the  development  and  implementation  of  a  special  nonlinear  solution  procedure 
consisting  of  two  steps: 

•  Design  load:  First  the  design  load  is  applied  to  the  bonded  joint  and,  after  a  material  nonlinear 
analysis  is  performed,  the  margins  of  safety  are  computed  for  the  different  modes  of  failure  identi¬ 
fied  for  the  joint  as  described  in  the  previous  section.  For  each  mode  of  failure  a  different  expres¬ 
sion  for  the  margin  of  safety  is  required,  which  may  be  based  on  average  quantities  as  described 
above. 

•  Failure  load:  After  the  design  load  analysis,  it  is  necessary  to  determine  the  failure  load  based  on 
the  margins  of  safety.  The  evaluation  of  failure  load  by  an  automatic  procedure  from  the  handbook 
framework  of  Stress  Check  required  developing  a  special  nonlinear  solution  procedure  which  per¬ 
forms  automatic  nonlinear  runs  with  load  increments  specified  by  the  user  and  checks  at  the  end  of 
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Integral  Average  Extraction 


(0)  Solution  *  SOL,  rune  #1  to  #8 

Total  length  of  integration  path  *  4.500000e-003 
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FIGURE  45.  Extraction  of  the  average  oy  along  an  element  edge. 

each  load  step  to  determine  whether  any  of  the  failure  criteria  have  been  exceeded  (MS  <  0).  Once 
the  first  failure  mode  is  reached,  it  is  possible  to  stop  the  analysis  or  continue  to  increment  the  load 
until  all  margins  are  exceeded. 

The  implementation  of  the  automatic  procedure  was  developed  using  the  material  nonlinear  capability 
of  Stress  Check  based  on  the  deformation  theory  of  plasticity  and  the  von  Mises  yield  criterion.  Only 
the  adhesive  material  have  elastic-plastic  properties  which  can  be  defined  as  elastic -ideally-plastic, 
bilinear,  Ramberg-Osgood  or  by  five  parameters  that  characterize  the  stress-strain  curve. 

Example:  Consider  the  bonded  lap  joint  shown  in  Figure  46.  Only  one  quarter  of  the  joint  is  consid¬ 
ered  for  the  2D  plane-strain  analysis  as  shown  in  the  figure.  The  first  few  layers  of  the  laminate  are 
included  explicitly  and  the  rest  are  lumped  into  a  sublaminate.  The  following  material  properties  are 
specified  for  the  composite  layers  and  the  adhesive: 
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T7“ 


P 


FIGURE  46.  Bonded  lap  joint. 


Composite:  £L=23.5xl06  psi,  £t=1.3x106  psi,  Gi7=0.55xl06  psi,  G^O^xlO6  psi,  vL7=0.3, 
Vtt=0.46 

Adhesive:  £=251,900  psi,  v=0.34,  gy=1 1,600  psi  (elastic-plastic). 

A  design  load  is  applied  to  the  bonded  joint  (P=2,000  lb),  and  a  material  nonlinear  analysis  is  per¬ 
formed  to  determine  whether  any  of  the  following  margin  criteria  are  exceeded  at  the  locations  indi¬ 
cated  in  Figure  47: 


MS,  =  l 

*  c  _ 


(67) 


Note  that  a  total  of  10  margin  criteria  are  considered  for  the  composite  joint,  5  on  the  left  side  and  5 
on  the  right  side.  The  failure  types  are  delamination  (Delam),  interface  failure  (Interf)  and  adhesive 
failure  (Adhes). 
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Delam 

Interf 

Adhes 

Extraction 

Min/Max 

Edge 

Edge 

Margin  of  Safety 

ms2 

ms2 

MSj 

Average 

NO 

YES 

YES 

Locations 

4 

4 

2 

FIGURE  47.  Extraction  location  and  type  for  limit  load  analysis. 

The  load  is  increased  (by  100  lb  increments)  and  the  nonlinear  analysis  continues  until  convergence  is 
achieved.  At  the  end  of  each  converged  solution  all  margin  criteria  are  checked  and  if  any  margin  is 
zero  or  negative  (MS  <  0),  the  magnitude  of  the  load  and  the  corresponding  margin  value  is  recorded. 
A  solution  record  is  also  kept  for  further  post-processing.  Figure  48  shows  the  settings  of  the  limit 
load  interface  in  Stress  Check. 

After  the  sequence  of  linear  solutions  is  available,  one  solution  is  selected  to  start  the  nonlinear  analy¬ 
sis.  The  type  of  nonlinear  run  (material  or  general)  is  selected  together  with  the  convergence  criteria 
(stress  or  energy),  allowable  tolerance  (%)  and  the  iteration  limit  for  each  load  increment.  Next,  the 
load  parameter  information  is  completed,  including  the  load  increment  and  the  upper  bound  for  the 
load.  It  is  possible  to  stop  the  analysis  after  the  first  margin  is  exceeded  by  checking  the  box  labeled 
‘stop  at  1st  M.S.  failure’.  If  left  unchecked,  the  solution  will  proceed  until  all  margin  criteria  are 
exceeded  or  the  upper  bound  of  the  load  is  reached,  whatever  occurs  first.  Finally,  all  the  relevant 
margin  criteria  defined  for  the  problem  are  selected  in  the  scroll  window  and  the  limit  load  analysis  is 
initiated. 
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Selected  margin 
criteria 


StressCheck  Limit  Load  Solver 


Linear  )  Nonlinear  ]  Model  |  Measurement  Limit  Load  |  SOLVE  1 1 
Solution  Run  Type  DOF  i:V"  '■  IrTTTT^  "■ 

A  .Lin.  ,  3008  3  . -1.  ~ 

,5  ,Lin.  ,  4484  Convergence:  [Energy  " 

.6  ,Un.  ,  6312 


SOLI 

SOLI 

SOLI 

SOLI 


,7  ,Lin.  ,  8492 


Lin.  .  11024 


Linear  Sol.:  SOLI 


F V  ft 

i  _ 

.  ,  hr 


Tolerance  (%):  j.25 
Iteration  Limit  |40 
Stop  at  1st  MS.  failure 


flH  nplaml  1!  Hi 

LJ  1  Cl  f  f  1  1 L.  mmm  B 

KB  DelamDL  HI 

r*  i  i  r* 

Param.:  p  tJ  Initial: 

2.000e^003  ^ 

L'Hlizt.mL'r.  Wg£; 

Step:  |l  00  j  Limit 

5000  ^ 

■■  UelamL’h  H| 

> 


Nonlinear 

settings 


Initial  load 

Upper  bound 
-  Load  increment 


FIGURE  48.  Limit  load  interface.  Setting  for  example  problem. 


After  the  limit  load  analysis  is  completed,  a  report  window  provides  a  summary  of  the  results  as 
shown  in  Figure  49.  Two  margin  criteria  were  exceeded:  An  interface  failure  occurred  on  the  left  side 


Plate 

Limit  Load  Summary 
with  central  hole.  NonL inear  model 

ID=  SOLI,  run  #8 

Criteria  Number 

Margin  Criteria 

Limit  Load 

Margin  Value 

1 

DELAMUL 

5 .000000e+003 

5.114385e-001 

2 

DELAMDL 

5 .000000e+003 

1.555231e-001 

3 

DELAMUR 

5 . 000000e+003 

9.346218e-001 

4 

DELAMDR 

5 .000000e+003 

6.299574e-001 

5 

ADHESL 

4 .000000e+003 

-2.606208e-002 

6 

AD HE SR 

5 .000000e+003 

8.395282e-001 

7 

INTERFUR 

5.000000e+003 

7.228068e-001 

8 

INTERFDR 

5 .000000e+003 

6.769590e-001 

9 

ZHTERFDL 

2.400000e+003 

-1.849117e-002 

10 

INTERFUL 

5 .000000e+003 

1.762295e-001 

2  out  of  10  Margin  criteria  exceeded. 

FIGURE  49.  Results  from  the  limit  load  analysis  of  the  example  problem. 


of  the  joint  just  below  the  adhesive  line  (InterfDL)  at  a  load  of  2,400  lb,  and  an  adhesive  failure  devel¬ 
oped  on  the  left  side  (AdhesL)  at  4,000  lb. 
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Finally,  Figure  50  shows  the  elastic-plastic  equivalent  (von  Mises)  stress  distribution  in  the  adhesive 
layer  for  P=2,400  lb,  P=4,000  lb  and  P=5,000  lb. 


P  =  2,400  lb 


P  =  4,000  lb 


P  =  5,000  lb 


FIGURE  50.  Elastic-plastic  equivalent  (von  Mises)  stress  distribution  following 
a  limit  load  analysis. 


StressCheck  V5J.afi 

NLMAT  ID  -  INTERFDL  5 
Run  *8.  OOF-1 1124 
Fnc.  •  Seq 
Max  ■  1.16006*004 
Min  ■  B. 6694e+001 


1.16009+004 
1.0771e+004 
9.9429e+003 
9.1 143e+003 
8.2857e+003 
7.457 1  e+003 
6.6286e+003 
5.8000e+003 
4.97146+003 
4.1429e+003 
3.3143e+003 
2.4857e+003 
1.8571e+003 
8.2857e+002 
0.0000e+000 


4.9  Task  9:  Sub-laminate  property  input 

A  capability  was  implemented  to  facilitate  the  input  of  orthotropic  material  properties  for  individ¬ 
ual  plies  and  for  sub-laminate  properties  for  2D  plane-strain  and  3D  analyses. 

•  Individual  plies :  When  the  ply-angles  are  not  contained  in  the  standard  2D  working  plane  (the 
plane-strain  XY  plane),  the  material  matrix  needed  for  the  2D  analysis  are  extracted  from  the  3D 
material  matrix.  The  material  coefficients  are  entered  in  the  material  directions,  and  then  assigned 
to  the  elements  together  with  the  ply  angle  information.  The  program  performs  all  necessary  trans¬ 
formations  to  compute  the  equivalent  2D  properties  in  the  XY  work  plane.  The  appropriate  3D 
properties  are  be  applied  when  the  planar  model  is  extruded  into  3D  as  well. 

•  Sub-laminate s:  When  a  set  of  plies  need  to  be  combined  in  a  single  layer  (sub-laminate),  the  prop¬ 
erties  of  the  sub-laminate  are  obtained  by  homogenization.  Again,  the  3D  material  coefficients  in 
the  material  axes  of  the  composite  are  entered  together  with  the  stacking  sequence  of  the  sub-lami¬ 
nate,  and  the  program  computes  the  equivalent  2D  properties  in  the  Stress  Check  XY  work  plane 
for  the  stack.  The  appropriate  3D  properties  are  applied  when  the  planar  model  is  extruded  into  3D 
as  well. 
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The  procedures  implemented  in  Stress  Check  to  account  for  these  two  activities  include  the  following: 

1  The  material  coefficients  of  a  single  orthotropic  ply  are  entered  in  the  direction  of  the  material  axes 

(Figure  51a).  Nine  engineering  coefficients  ( E\\ ,  f?22>  ^33»  ^23*  ^3i>  V12»  v23>  vn)  the 

direction  of  the  material  axes  and  the  ply  thickness  must  be  provided.  The  three  coefficients  of 
thermal  expansion  (ay,  (X22,  0133)  and  the  mass  density  can  also  be  entered,  but  they  are  not 
required  unless  there  is  thermal  loading  or  a  modal  analysis  is  needed. 

2  After  the  material  properties  for  a  single  ply  are  entered,  the  ply  group  information  must  be  pro¬ 
vided.  This  includes  the  angular  orientation  of  each  ply  in  the  group,  in  accordance  with  the  fol¬ 
lowing  convention:  A  positive  ply  angle  (0)  is  measured  as  a  counterclockwise  rotation  about  the 
z-axis  of  a  local  coordinate  system  (xyz).  The  z-axis  of  the  local  is  aligned  with  the  material  3-axis 
as  shown  in  Figure  51b. 


(a)  Material  axes 


(b)  Rotation  of  material  axes 


FIGURE  51.  Material  coordinate  systems  for  orthotropic  plies. 


3  The  ply  groups  are  assembled  into  a  stack  by  providing  the  total  number  of  layers  in  the  laminate, 
defining  whether  the  stack  is  symmetric  or  not  and  entering  the  stacking  sequence  based  on  the  ply 
group  names.  This  procedure  provides  a  great  flexibility,  since  various  ply  groups  of  different 
material  properties  can  be  assembled  together  to  represent  a  single  sub-lamina. 

4  Finally,  the  stack  is  assigned  to  the  elements  in  the  mesh  identifying  the  local  coordinate  system, 
the  z-axis  of  which  is  assumed  to  coincide  with  the  material  3-axis  and  must  be  perpendicular  to 
the  plane  of  the  ply  (Figure  51b).  The  case  of  a  single  ply  assigned  to  an  element  is  treated  as  a  par¬ 
ticular  case  of  a  sublaminate  with  one  orthotropic  layer. 

The  two-dimensional  material  stiffness  matrix  needed  for  the  plane-strain  analysis  is  then  extracted 
form  the  three-dimensional  material  matrix  in  global  coordinates  (XYZ).  When  a  group  of  laminae  is 
assigned  to  a  single  element,  the  three-dimensional  effective  properties  are  computed  by  homogeniza¬ 
tion,  following  the  procedure  described  in  Ref.  [19]  and  discussed  in  detail  below. 
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4.9.1  Formulation 

The  3D  stress-strain  relations  for  an  orthotropic  material  in  the  material  coordinate  system  (Figure 
51a)  can  be  written  as: 

{o}  =  [C]{e>  (68) 

where. 


C\ i  C ]2  C13  0  0  0 

C22  C23  0  0  0 

[C]  =  C33  0  0  0 

sym.  C#  0  0 

C44  0 
£55 


(69) 


[C]  is  the  material  stiffness  matrix  with  nine  independent  constants.  The  relation  between  the  Cy  and 
the  engineering  coefficients  is  given  by  (Ref.  [20],  page  41): 


C]1  —  (1  -  ',23V32)V£ji, 

C,2  =  (v21  +  v23v3l)^^ll’ 


C 22  ~  (  1  ~  V31  v13)^^22» 

C,3  =  (v,3  +  v23v12)V£33, 


C33  -  (1 -v12v21)V£33 

^23  =  (v32  +  V12V3l)^22 


£44 

r*“i 

CM 

II 

^5?  =  &319 

^66  ~  @12 

V  = 

0-Vl2V2l 

~  V23V32“  V31V13" 

-2v12v23v31)_1 

E 22 

^33 

E33 

=  V 

12  r  ’ 

E\l 

V32  “  V23^“  * 
*■'22 

V31  ~  V13p 
^■11 

Next  consider  the  situation  shown  in  Figure  5 1  b,  in  which  a  rotation  0  about  the  z-axis  (coincident 
with  the  material  3-axis)  is  performed  for  the  fcth  layer  in  the  laminate.  The  stress-strain  relations  in 
the  xyz  local  system  can  be  written  as: 

{°hyz  =  lQ(k)UehyZ  (70) 

where. 


54  of  68 


Stress/Failure  Analysis  Software  for  Multi-Material  Interfaces 


Implementation  and  results 


Wxvz  = 


{e}„z  = 


lQ(k)]  = 


Q\\  Q12  Qn  Q\(>  o  o 

622  ^23  G26  0  0 

Q33  Q 36  0  0 

sym.  0  0 

Q 44  645 

e55 


The  material  matrix  for  the  fcth  layer  in  the  local  system  [Q(k)]  is  obtained  by  the  following  transforma¬ 
tion  (Ref.  [20],  page  31): 

[Q(k)]  =  [Tf[C][T] 


m  n  0  0  0  mn 

n  m  0  0  0  -mn 

0  0  1  0  0  0 

0  0  0  m  -n  0 

0  0  0  n  m  0 

-2 mn  2 mn  0  0  0  m  —n 


where  [C]  is  given  in  Eq.  (69),  m  =  cos0  and  n  =  sin0.  Equations  (70)  and  (71)  represent  the  stress- 
strain  relations  for  a  lamina  in  a  local  coordinate  system  after  a  rotation  0  about  the  z-axis. 

Next  we  consider  a  laminate  of  constant  thickness  h  composed  of  AMayers  of  thickness  t \ 
(h  =  S'"  tk).  Each  lamina  can  have  an  arbitrary  orientation  (z-rotation)  with  respect  to  the  local 
coordmate'system  (xyz).  The  properties  of  this  laminate  can  be  effectively  represented  by  an  homoge¬ 
neous  anisotropic  solid.  Consider  the  stress-strain  relations  for  the  homogeneous  solid  in  the  local 
system  of  the  form: 


1 Q )  = 


1°},VZ  = 

IfiKe} 

xyz 

fi.1  fi>2 

fil3 

fi.6 

0 

0 

Qn 

Ql3 

Q2 6 

0 

0 

Q33 

Q3 6 

0 

0 

(72) 

sym. 

Q(& 

0 

0 

^44  ^45 

e5i 


Stress/Failure  Analysis  Software  for  Multi-Material  Interfaces 


55  of  68 


Implementation  and  results 


where  the  effective  elastic  constants  of  the  solid  are  computed  from  the  material  matrices  of  the 
individual  plies  by  the  transformations  (Ref.  [19]): 


Once  the  material  matrix  for  the  homogenized  solid  is  computed  in  the  local  system  (xyz),  we  need  to 
determine  the  material  matrix  in  the  global  system  (XYZ)  shown  in  Figure  51b,  and  then  extract  the 
material  matrix  associated  with  the  ‘plane-strain’  plane  (XF-plane). 
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In  the  global  system,  the  three-dimensional  stress-strain  relations  can  be  written  in  the  following 
form: 

{°}xrz  =  [^3DHe}xrz  (7*) 

where  the  global  material  matrix  [E^p]  of  the  homogenized  solid  is  computed  from  the  local  matrix 
[Q]  as  follows: 

[Ew]  =  [T0]t[Q][T0]  (75) 

The  transformation  matrix  [Tq]  is  given  by: 


1 T0]  = 


£ 

£ 

v2 

hh 

/3/i 

2 

2 

2 

*1 

m2 

m3 

mxm2 

m2m3 

m3mx 

2 

2 

2n 

»i 

n2 

n\n2 

n2n3 

n3rtl 

2/|  YYl  | 

2  l2m2 

2/3/W3 

l\m2  +  hm\ 

/2/W3  *4"  l3ftl2 

/3mj  +  /j/n3 

2mxnx 

2  m2n2 

2/n3n3 

mxn2  +  m2nx 

m2n3  +  m3n2 

^3^!  +  m  x  n3 

2nxlx 

2/22^2 

2  W3/3 

nxl2  +  n2lx 

n2l3  "4"  W3/2 

»3/l+»l/3. 

where  m,-,  nt-  are  the  direction  cosines  that  relate  a  point  in  the  global  ( XYZ)  coordinate  system  to  a 
point  in  the  local  (ryz)  coordinate  systems: 


(i  h  h 

m  i  m,  m, 


Since  the  transformation  indicated  by  Eq.  (75)  is  very  general,  the  3D  global  material  matrix  can  be 
fully  populated.  Therefore,  the  material  matrix  in  Eq.  (75)  can  be  written  as  (using  a  single  index 
notation  for  the  21  independent  coefficients): 


[E3D]  = 


E\  E2  Ea  Ef  £u  E16 

£3  E5  £g  En  £]7 

E6  £9  £, 3  £16 

sym.  £)0  £14  £,9 
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The  2D  plane-strain  material  matrix  is  extracted  for  Eq.  (78)  by  pulling  out  the  rows  and  columns 
associated  with  the  XT-plane  of  the  global  system: 

ax  Et  E2  E7  ex 

oY  =  £3  £g  ey  <79> 

[sym  £ioJ  LYxr 

The  material  matrix  in  Eq.  (79)  is  used  for  the  computation  of  the  stiffness  matrices  for  2D-elasticity 
under  plane-strain  conditions,  and  the  material  matrix  in  Eq.  (78)  is  used  for  the  computation  of  the 
stiffness  matrices  for  3D-elasticity  for  those  elements  with  laminate  property  assignments. 

4.9.2  Examples 

Two  examples  are  presented  to  illustrate  the  implementation.  In  the  first  example  each  layer  of  the 
laminate  is  explicitly  included  in  the  model,  while  in  the  second  example  some  plies  are  replaced  by  a 
single  element  with  homogenized  properties.  A  comparison  between  the  two  is  performed  to  assess 
the  influence  of  the  homogenization  in  the  results. 

Example  1:  Consider  a  16-ply  [4(0/90)]$  simply  supported  laminated  composite  strip  under  plane 
strain  conditions  subjected  to  a  sinusoidal  traction  as  shown  in  Figure  52.  Each  ply  is  of  the  same 


FIGURE  52.  16  ply  [4(0/90)]s  simply-supported  laminated  composite  under  sinusoidal  load.  Notation  and  mesh. 

material  and  thickness  (r=0.0625  mm)  with  the  following  properties: 


58  of  68 


Stress/Failure  Analysis  Software  for  Multi-Material  Interfaces 


Implementation  and  results 


El  =  1.38xl05  MPa,  ET  =  9.3xl05  MPa 
Glt  =  4.6x10s  MPa,  Gjj  =  3.1x10s  MPa 
v^7’=  0.3,  Vjj—  0.5 

Because  of  symmetry,  only  half  of  the  beam  was  considered  for  the  analysis.  The  mesh  shown  in  Fig¬ 
ure  52  consists  of  16  quadrilateral  elements,  one  for  each  ply.  The  material  properties  for  a  each  ply 
are  entered  in  the  material  definition  form  shown  in  Figure  53a.  All  nine  engineering  coefficients 
(*n  ,  E2 2,  £33,  Gj 2,  G2 3,  G3],  Vj2,  V23,  v j 3)  in  the  direction  of  the  material  axis  and  the  ply  thickness 
must  be  provided  in  this  form.  The  three  coefficients  of  thermal  expansion  (a^,  a22,  (X33)  and  the 
mass  density  are  optional.  After  the  material  properties  for  a  single  ply  are  entered,  the  ply  group 
information  must  be  provided  as  shown  in  Figure  53b:  The  ply  group  name,  the  lamina  material  name 
and  the  ply  layout  angles.  The  stack  form  shown  in  Figure  53c  is  used  to  assemble  the  ply  groups  by 
providing  the  stack  sequence  name,  the  number  of  layers  in  the  laminate,  whether  the  stack  is  sym¬ 
metric  or  not  and  the  stacking  sequence  based  on  the  ply  group  names.  As  shown  in  Figure  53,  the 


Math  j  Sadi  on  Prop.)  ThicknaM  Materia]  j  u»<  |  >  P^  |  Slack  | 

iDefine  HlUnaar  *»1|  Selection  "»1  Ply  Name  Ptopeiy  Loyo  it  Angles 


Material  ID  Type 

|SINGli^'^™,ljminate_Ortt)0 
T  ID:  jsiNGLEPLY  : 
r  Scale:  1 0.000e*000 


~j  X  Rt 
its:  (US 


Option:  (Defined  Mtrl.  Units:  (US  Zj 

Material  | Linear  H  Fitlmfl:  |no 

Type:  (Lamin  Ortho  "z\  C®8*  |  R  Stress 
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j1.3B0e+005 
|9.300e+003 
|5.000e-001 
|4  600e+003 
|4  600e*003 
|0.000e*000 


E2  2:  |9.300e+0Q3 
v!2:  [3  OOOe-OoT 

v13:  ]3  OOOe-OOl 
G23:  |3.100e*003 

Dens.:  |l  000e-004 
822:  (O.OOOe+OOO 


Replace 

Delete 

Cancel 

* 

(a)  Material  definition 


Ply  Group  Name:  |R_Y90 

Umina  Materia)  Name:  JsiNGLEPLY 
Ply  layout  angles:  (e.g.  (M5H530) 

I5  — — 

_ leave  blank  » ISOTROPIC 

I  Accept  I  Replace!  Delate  I  Pui 


(b)  Ply  group  definition 


(c)  Stacking  definition 


Ply  Stack  | 
Stack  Name 
Ilaho 


Layers  Sequence 

*7  1,  plyO 


Stack  Sequence  Name:  |LAM90 
Number  of  Layers:  |i  )Non-Symmetric^J 
Stacking  sequence:  (e.g.  ocb*3:2*c) 

jpiyoo 

_  (e  g.  a  b,  c  ore  piy  group  names) _ 

1  Accept  )  Replace  I  Delete  I  Purge  J 


FIGURE  53.  Input  forms  for  laminate  properties  definition.  Example  1 
material  property  SINGLEPLY,  defined  in  the  material  definition  form,  is  used  to  define  the  ply  group 
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PLY90  with  a  ply  layout  of  90  degrees  with  respect  to  the  z-axis  of  a  local  system  and  the  stack 
LAM90  is  created  with  a  single  ply  using  the  plygroup  PLY90.  The  last  step  is  to  assign  the  material 
stacks  to  the  elements  and  select  the  local  coordinate  system,  the  z-axis  of  which  is  assumed  to  be 
aligned  with  the  material  3-axis.  The  material  assignment  form  is  shown  in  Figure  53.  The  orientation 
of  the  local  coordinate  system  relative  to  the  global  system  is  used  to  determine  the  material  matrix  in 
the  global  system  as  explained  in  the  previous  section. 


Stack  name 


Local  system 


FIGURE  54.  Input  form  for  laminate  properties  assignment  to  elements. 


A  sequence  of  finite  element  solutions  was  obtained  for  polynomial  orders  ranging  from  1  to  8  for  a 
length  to  thickness  ratio  ( a/h )  of  20.  The  estimated  relative  error  in  energy  norm  and  the  deformed 
configuration  for  run  #8  are  shown  in  Figure  55.  The  largest  vertical  deflection  of  the  laminated  strip 
occurs  along  the  symmetry  line  (Uy  =  -0.24125). 

The  distribution  of  the  normal  stress  ax  is  shown  in  Figure  56  for  run  #8.  The  maximum  value  occurs 
at  the  lower  surface  along  the  symmetry  line  (cxmax  =  394).  These  results  will  be  compared  with  those 
of  example  2. 
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Error  Estimate 
16-ply  lminate  [4(0/90)]s 
(0)  Solution  -  SOL,  runs  Ml  to  MB 


Run  M 

DOF 

Total  Potential  Energy 

Convergence  Rate 

%  Error 

1 

34 

-3 . 6773 39739 17 le-001 

0.00 

93.70 

2 

100 

-2 . 72  6732225 119e+000 

1.03 

30.75 

3 

166 

-3 . 004822 685327e+000 

3.69 

4.73 

4 

2  64 

-3 . 0 1 1498754148e+000 

4.84 

0.50 

5 

394 

-3 . Q11573787861e+000 

6.25 

0.04 

6 

556 

-3 . 011574291076e+000 

7.83 

0.00 

7 

750 

-3 . 011574292972e+000 

2.92 

0.00 

8 

976 

-3 . 011574293285e+000 

2.92 

0.00 

Estimated  Limit  -3 .0115742933 69909e+000 


StressCheck  VS.O.bM 

LINEAR  ID  •  SOL 
Run  -  6  ,  DOF-976 
Deformed 

Hex  -  1. 8065e-002 

Min  -  -2.4125e-<X)l 


FIGURE  55.  Estimated  relative  error  in  energy  norm  and  deformed  configuration  for  example  1. 


StressCheck  V5.0.bl4 

LINEAR  ID  -  SOL 
Run  -  6  ,  DOF-976 
Fnc.  -  Sx 

Hex  -  3.9422^002 

Hin  -  -3.9422^002 

3.9422et002 
3.1538e*002 
2.3653^002 
1.5769fr+002 
7.8844e+001 
-1.5259e-005 
-7.8844C+001 
-1. 5769e*002 
-2.3653e-t002 
-3.1538e+002 
-3.9422frtOQ2 


FIGURE  56.  Normal  stress  distribution  (Sx)  for  example  1. 
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Example  2:  Consider  the  same  16-ply  [4(0/90)]5  simply  supported  laminated  composite  strip  of 
example  1,  but  with  the  eight  central  layers  replaced  by  a  single  element  with  equivalent  properties. 
The  finite  element  mesh  is  shown  in  Figure  57.  The  ply  group  definition  and  stacking  sequence  for  the 


FIGURE  57.  Mesh  for  a  16  ply  [4(0/90)]s  simply-supported  laminated  composite  under  sinusoidal 
load  with  the  eight  central  layers  replaced  by  a  single  element. 

8  individual  layers  (4  above  and  4  below  the  sub-laminate)  are  the  same  as  shown  in  Figure  53  for 
example  1.  The  corresponding  input  records  for  the  sub-laminate  are  shown  in  Figure  58.  The  central 
8-plies  [2(0/90]s  can  be  described  in  more  than  one  way  in  the  input  field  of  Figure  58a.  For  example. 


Laminated  Material  Definition 


Pfy  |  Stack  |  ' 

Ply  Mama  :  Propary  Layout  Anglur 


Ply  Group  Noma:  | 

Lamina  MatarialNama:  J  LAMINA 

layout  ongiaa:  (i.g.  E4S>45;3flj) 


2[0/90];2[90/0] 


Laovt  Wank  If  ISOTROPIC 


Accapt  pRaplaca  |  Dalata  |  Purga 


Laminated  Material  Definition 


PV  Stack 
Stack  Nama 


i  ijayarft  Sequence 


Stack  Saguanca  Nama:  jSUBLAMINA 
Numbaroflayam:  je  :|Nort-Symmetricj^j 
Stacking  saguenos:(a,g.  a:b*3;2*c) 


F 


sublam 


(e.g.gix  caia  ply  group  namas) 


Accept  1  ftepiaca|  Delate  |  Purga 


(a)  Ply  group  definition 


(b)  Stacking  definition 


FIGURE  58.  Input  forms  for  sub-laminate  properties  definition.  Example  2 


0/90/0/90/90/0/90/0,  or  2[0/90];2[90/0]  as  shown.  Alternatively,  define  only  half  of  the  sub-laminate 
(4-plies)  in  the  ply  definition  form  (that  is,  2[0/90])  and  then  use  the  Symmetric  option  in  the  stacking 
sequence. 

A  sequence  of  finite  element  solutions  was  obtained  for  polynomial  orders  ranging  from  1  to  8  for  a 
length  to  thickness  ratio  ( a/h )  of  20.  The  estimated  relative  error  in  energy  norm  and  the  deformed 
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configuration  for  run  #8  are  shown  in  Figure  59.  Note  that  the  values  of  the  total  potential  energy  for 
run  #8  (-3.1 12472)  and  largest  vertical  deflection  for  run  #8  (-2.4933)  are  very  close  to  those  of  exam¬ 
ple  1  (-3.011574  and  -2.4125,  respectively). 

Error  Estimate 
16-ply  lminate  [4(0/90)]s 
(0)  Solution  -  SOL,  runs  #1  to  #6 


FIGURE  59.  Estimated  relative  error  in  energy  norm  and  deformed  configuration  for  example  2. 


The  normal  stress  distribution  ax  for  run  #8  is  shown  in  Figure  60.  Note  again  that  the  maximum 


stress  occurs  at  the  lower  surface  along  the  symmetry  line  (axmax  =  408)  and  it  is  very  close  to  the 
value  obtained  for  example  1  (axmax  =  394). 
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The  analysis  was  repeated  for  several  length-to-thickness  ratios  (a/h)  and  the  results  are  summarized 
in  Table  2.  As  can  be  seen,  the  results  are  very  close,  demonstrating  the  good  approximating  charac- 


TABLE  2.  Results  for  examples  1  and  2 


a/h 

Example 

Potential 

Energy 

»t  min 
uy 

<r  max 

ux 

20 

394 

408 

40 

1 

-92.6603 

BH 

1569 

2 

-95.9683 

mBm 

1626 

80 

1 

-2935.45 

-58.714 

6267 

2 

6495 

teristics  of  the  homogenization  procedure  implemented  for  the  sub-laminate. 

Both  example  problems  were  analyzed  again,  but  this  time  in  a  fully  three-dimensional  setting  using 
the  extrusion  option  in  Stress  Check.  Because  the  composite  strip  is  a  cross  ply  laminate,  the  results 
should  be  identical  to  those  of  plane-strain,  provided  that  the  normal  displacements  of  the  faces  con¬ 
tained  in  the  XY-plane  are  constrained. 

In  both  cases  a  sequence  of  finite  element  solutions  was  obtained  for  polynomial  orders  ranging  form 
1  to  8.  Figure  61  shows  the  estimated  relative  error  in  energy  norm  for  both  cases.  Note  that  the  poten¬ 
tial  energy  values  are  the  same  as  for  the  plane-strain  cases  (see  Figures  55  and  59),  but  the  number  of 
degrees  of  freedom  (DOF)  have  increased  substantially  for  the  3D-solution. 


GlobalError  GlobalError 

16-ply  Iminate  [4 (0/90)] s  16-ply  laminate  [4 (0/90)] a 

(0)  Solution  -  SOL/  runs  #1  to  00  (0)  Solution  *  SOL,  runs  #1  to  #8 


Run  * 

DOF 

Total  Potential  Energy 

Convergence  Rate 

%  Error 

Run  0 

oor 

Total  Potential  Energy  Convergence  Rate 

%  Error 

1 

68 

-3 . 6773 3973 9127e-001 

0.00 

93.70 

1 

40 

-3 .754066833841e-001 

0.00 

93.78 

2 

251 

-2.726732224908e+000 

0.85 

30.7$ 

2 

146 

-2 . 81235916882 6C+000 

0.85 

31.05 

3 

434 

-3 . 00462  2  685147e+000 

3.42 

4.73 

3 

252 

-3 . 105194949440e+000 

3.41 

4.84 

4 

760 

-3 . 0114987 67466e+000 

3.83 

0.50 

4 

451 

-3 . 1123789362 lle+000 

3.75 

0.55 

5 

1269 

-3.011573787270e+000 

5.11 

0.04 

5 

743 

-3 . 11247079141164000 

4.96 

0.05 

6 

2009 

-3 . 011574290980e+000 

2.25 

0.01 

6 

1155 

-3 . 112471534477e+000 

2.27 

0.02 

7 

2988 

-3 . 0115742 93356e+000 

0.00 

0.01 

7 

1714 

-3. 11247153863 le4O00 

0.00 

0.02 

6 

4274 

-3 . 0115742 92866e4O00 

0.01 

0.01 

8 

2447 

-3 . 112 47153 85 65e4000 

0.00 

0.02 

Estimated  Limit  -3 ,011574230932173e+000  Estimated  Limit  -3 . 112471446182672e+000 


(a)  Example  1  (b)  Example  2 


FIGURE  61.  Estimated  relative  error  in  energy  norm  from  the  3D  analysis  of  examples  1  and  2. 
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Figure  62  shows  the  deformed  configuration  for  mn  #  8  for  both  cases.  Again  the  largest  displacement 
is  the  same  as  in  the  case  of  plane-strain  (compare  with  Figures  55  and  59). 


FIGURE  62.  Deformed  configuration  from  3D  analysis  of  examples  1  and  2. 


Finally,  Figure  63  shows  the  normal  stress  distribution  ox  for  run  #8  for  both  examples,  which  are 
identical  to  their  plane-strain  counterpart  (compare  with  Figures  59  and  60). 


FIGURE  63.  Normal  stress  distribution  (Sx)  from  the  3D  analysis  of  examples  1  and  2. 
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4.10  Task  10:  Delivery  of  Stress  Check  software 

One  copy  of  Stress  Check  (Release  5.0)  software  executable  on  a  Windows  NT  workstation  and 
pertinent  documentation  will  be  delivered  to  Dr.  Douglas  J.  Holzhauer,  Mechanical  Engineer,  Air 
Force  Rome  Laboratory,  RL/ERDS,  Griffiss  AFB,  NY  13441-5700.  Rome  Laboratory  shall  have  the 
same  rights  to  Stress  Check  as  a  “paid-up”  licensee. 

All  technical  capabilities  described  in  this  report  were  incorporated  in  Release  5.0  of  Stress  Check 
which  is  currently  undergoing  beta  testing.  The  official  release  date  is  December  1st  1999,  at  which 
time  the  software  will  be  delivered  to  the  Air  Force. 

5  Summary  and  conclusions _ 


The  technological  developments  implemented  during  the  Phase  II  project  are  an  essential  prereq¬ 
uisite  to  proper  interpretation  of  experimental  data  in  much  the  same  way  as  the  ability  to  compute 
stress  intensity  factors  in  linear  elastic  fracture  mechanics  is  an  essential  prerequisite  to  proper  inter¬ 
pretation  of  experimentally  obtained  crack  propagation  data.  The  new  capabilities,  coupled  with  refer¬ 
ence  data  obtained  from  simple  experiments,  will  make  it  possible  to  evaluate  design  alternatives  in 
the  fields  of  electronic  component  design  and  composite  materials  technology.  This  capability,  made 
available  for  professional  use  through  the  finite  element  analysis  software  (Stress  Check),  is  expected 
to  be  of  substantial  interest  to  manufacturers  of  electronic  components,  aerospace  companies,  and 
their  suppliers. 

A  change  in  the  research  objectives  was  needed  because  during  the  course  of  our  investigation  an 
improved  and  more  generally  applicable  methodology  for  failure  analysis  of  composite  materials  was 
conceived.  The  new  technological  capability  provides  improved  numerical  simulation  tools  for  the 
design  and  analysis  of  laminated  composite  materials.  The  essential  differences  are  that  whereas  the 
original  task  description  envisioned  failure  analysis  of  composite  materials  entirely  by  means  of  lin¬ 
ear  models,  utilizing  the  concept  of  generalized  stress  intensity  factors,  the  new  tasks  advanced  the 
development  of  computational  support  for  failure  criteria  based  on  linear  as  well  as  nonlinear  material 
properties. 

The  modification  to  the  original  goals  of  the  Phase  II  project  addressed  the  requirements  of  a 
working  group  of  a  government-industry  consortium,  called  the  Composite  Affordability  Initiative 
(CAI).  This  group  was  charged  with  the  responsibility  of  identifying  state-of-the-art  numerical  simu¬ 
lation  technology  that  could  be  used  in  the  design  and  analysis  of  composite  joints  with  particular  ref¬ 
erence  to  the  development  of  next  generation  fighter  aircraft.  The  group  conducted  a  detailed 
evaluation  of  the  technology  we  have  developed  and  concluded  unanimously  that  Stress  Check  is  the 
best  available  tool  for  analyzing  and  designing  bonded  composite  joints.  The  evaluation  also  identi¬ 
fied  development  tasks  that  needed  to  be  completed  in  order  to  make  Stress  Check  fully  useful  for  the 
purposes  of  composite  joint  analysis.  The  new  technological  capabilities  implemented  as  described  in 
tasks  7  to  9  satisfies  those  additional  requirements.  They  provide  the  numerical  simulation  tools  that 
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can  be  used  in  the  design  and  analysis  of  laminated  composite  joints,  making  advanced  technology 
available  for  design  engineers  and  analyst  in  the  aerospace  industry. 
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